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ABSTRACT 


A computer code is developed and tested which simulates the transient 
evaporation of a single liquid droplet from the surface of a semi-infinite solid 
subject to radiant heat input from above. For relatively low temperature 
incident radiation, it is shown that the direct absorption of radiant energy by 
the droplet can be treated as purely boundary conditions, while a model for 
higher temperature incident radiation would require the addition of constant 
heat source terms. The heat equation is numerically coupled between the 
liquid and solid domains by using a predictor-corrector scheme. Three one- 
dimensional solution schemes are used within the droplet: a start-up semi- 
infinite medium solution, a tridiagonal Crank-Nicholson transient solution, 
and a steady-state solution. The solid surface temperatures at each time step 
are calculated through careful numerical integration of an axisymmetric 
Green’s functions solution equation with the forcing function given by the 
past lower droplet surface and solid-vapor boundary heat fluxes. The time 
step is increased after a sensitive initial period to allow for reasonable run 
times. Two geometry models are included which give the droplet height as a 
function of current droplet volume and initial wetted radius; the second 
allows inclusion of the effects of initial contact angle and receding angle. 
Using water as the liquid and Macor, a low-thermal conductivity material, as 
the solid, the program output was compared to the experimental results in 
this line of research. They correlate well to the experiments in which the 
critical geometric shape factor and evaporation time were most easily 


measured. 
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FOREWORD 


This report describes the research performed during the period September 1992 - August 
1993 under a joint research program between the Mechanical Engineering Department of 
the University of Maryland at College Park and the Building and Fire Research Laboratory 
of the National Institute of Standards and Technology. The research was conducted in the 
laboratories of the BFRL by G. White and S. Tinker, Graduate Research Assistants of the 
ME Department under the joint supervision of Dr. M. di Marzo (ME Department - UMCP) 
and Dr. D. Evans (BFRL - NIST). This report also constitutes the Master Thesis of Mr. G. 
White, which has been completed and will be defended in the month of December 1993. 
Ms. S. Thinker performed experiments on dropwise evaporation which will be included in 
the final report for next grant period. She was also responsible for the formulation of a two- 
parameters model for the description of the droplet shape which was used by G.White in 
his thesis. Ms. Tinker is currently working on the formulation and validation of a multi- 


droplet model for the prediction of sparse spray evaporative cooling. 
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1. INTRODUCTION 

Because of the very high latent heat of vaporization of water, cooling using — 
water droplets has many engineering applications. These uses vary from the 
quenching of metals and cooling of turbine blades to environmental control 
systems. Therefore, cooling of hot surfaces has been the subject of numerous 
investigations. Toda [1] and Bonacina [2] designed spray cooling experiments 
while more recent work by Inada [3] and Makino [4] studied the evaporation 


of single droplets heated from below by conduction. 


The theoretical transient single-droplet evaporative cooling computer code 
developed by the author was completed under a joint research program 
between the Mechanical Engineering Department of the University of 
Maryland and the Center for Fire Research of the National Institute of 
Standards and Technology. This research has focused on extending the 
current body of knowledge on evaporative cooling to the fire safety field. 
Previous experimental work [5-9] has investigated the temperature response 
of the surface of a semi-infinite solid surrounding a gently deposited water 
dropiet. The surface may be heated from below by conduction or from above 
by radiation, but previous theoretical work [8, 10-12] has focused on the 
conduction case. More recent work has focused on measuring [13] or 
predicting [14] the surface temperatures of a solid impacted by a random 
droplet array. The time for reaching the surface-averaged steady-state 
temperature for a surface subject to a sudden droplet array is critical to the 
protection of equipment exposed to radiant heat. Therefore, the results of the 
single-droplet code presented here will be used in a multi-droplet code. To 


reach that goal, the code documented here is designed to output transient 


solid surface temperatures and droplet evaporation times for water droplets 
of variable size and initial wetted area gently deposited on low-thermal 


conductivity semi-infinite solids subject to radiant energy. 


1.1 Problem Description 

The computer code documented here models the specific experimental 
geometry of Dawson [13] as shown in Figure 1. Heating is from above by 
three radiant panels operating in a temperature range of from 475-650°C. The 
droplet is gently deposited at the center of a semi-infinite solid Macor tile, a 
glass-like ceramic resistant to thermal stresses and chosen for its low thermal 
conductivity (1.29 W/m-K) and high emissivity of 0.84 [8]. The temperature 
profile along a line passing through the droplet center was measured using 
infrared thermography data recorded on video tape. Only initial surface 
temperatures less than that required for nucleate boiling in the droplet were 
allowed (i.e., Tso < 162°C). Therefore, the droplet vaporization process is 


exclusively evaporative. 


The droplet-solid geometry is shown in Figure 2. An overall energy balance 
predicts the amount of energy required to evaporate the droplet but cannot 
predict the relative proportions of the heat transfer mechanisms nor the 
evaporation time. Energy enters the droplet by conduction from below and 
radiation from above and leaves the droplet with the mass evaporated at the 
upper surface. A small amount of energy also leaves the upper surface by 
natural convection and is quantified using the experimentally measured 
value of h (approximately 10 W/m-K). The Chilton-Colburn analogy will be 
used to develop a relationship for the energy leaving by evaporation. 


Chapter 2 is devoted to quantifying the energy added to the droplet by 
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radiation. It is shown that the radiation may be assumed to be input 
exclusively at the droplet boundaries. The conductive contribution from 
below, responsible for the surface cooling effect, requires coupling of the heat 
equation in the liquid and solid domains. This is accomplished using a 
simple but effective predict-correct method. Surface-tension-induced motion 
in the droplet can be neglected because of the very small induced velocities, 
and Rayleigh-instability-induced flow is not of concern because of the large 
response times. Then, the thin film assumption that conductive heat transfer 
in the droplet is one-dimensional in the axial direction is made based on 
previous results [15] for the case of heating by conduction from below. This 
assumption is true everywhere and at all times except near the droplet edge 
where the radial flux component may be as large as 10% of the total flux. The 
radiation case droplet tends to be thinner than the conduction case one, so the 
assumption can safely be applied to the present problem. The solid domain 


is best handled by a Green's function Boundary Element Method (BEM). 


1.2 Key Assumptions 

The major assumptions used by the computer code are listed in the following: 

(1) Liquid flow and thus convection inside the droplet are negligible. 

(2) | Conduction in the droplet is one-dimensional in the axial direction. 

(3) A semi-infinite solid solution can sufficiently describe the early solid- 
liquid interfacial fluxes. 

(4) The direct radiation absorption by the droplet can be sufficiently 
described by radiation absorption at the boundaries only. 

(5) | The droplet shape may be empirically described by geometry Model A 
or B (subsection 3.9): spherical segment shape or pancake shape with 


possible shrinkage of the wetted area. 


2. SECTION I - RADIATIVE HEAT INPUT MODEL 

The evaporative process cools a surface by conducting sensible energy from 
the bulk of the solid in order to supply latent heat of vaporization. However, 
incident radiative heat flux (e.g. from a nearby fire) can also supply the 
needed energy and thus limit the cooling effect. The manner in which the 
radiative flux enters the system must affect the transient characteristics of the 
evaporative process. For example, radiation absorbed right at the droplet 
surface would affect the temperature-sensitive evaporative heat flux more 
strongly than radiation absorbed directly by the solid at the solid-liquid 
interface. Also radiation incident to the droplet at shallow angles tend to be 
reflected away from the solid entirely. Therefore, this section presents a 
radiative heat input model of enough detail so as to specify locations at which 
the radiation enters the droplet-solid system but not so complex as to make 
ti. problem extremely computationally intensive. This model is necessary 
for the program developed in Chapter 3, which simulates the controlled 
laboratory environment, because the overall energy balance provided by the 
known initial heat flux in the tile only gives the magnitude of the radiation 
absorbed by the tile but not the directions from which the radiation is 


incident. 


2.1 Derivation of Model Equation 

Because of the assumed azimuthal symmetry of the droplet and semi-infinite 
solid, the polar angle distribution of the incident radiation and not the 
azimuthal distribution is needed. Then, flux incident from polar angle @ with 
wavelength A can be designated E,,9. If one makes the assumption that the 


liquid-vapor interface is horizontal and flat, then the radiative flux is only a 


function of the depth z below the interface and not also of the radial position 
and current droplet shape. The droplet surface normal need not be 
calculated as a function of radial position and time. This assumption may be 
reasonable because the droplet is usually rather thin. Two additional 
assumptions are also made now: radiation scattering within the droplet is 
negligible and radiation reaching the liquid-solid interface is completely 
absorbed by the solid. Then with the effects of projected area, surface 
reflection, and absorption inside the droplet, the heat flux at a depth z below 


the droplet surface due to Ej 9 is 
F, 9(Z) = E, 9 cos@ (1 - pg) exp(-K,Z/[) (2.1) 


where ky is the spectral absorption coefficient of the liquid, p is the direction 
cosine to the z-axis, and Pg is the reflectivity of the liquid surface for an angle 
6 between the incident radiation and the droplet normal. x, acts to 
exponentially damp the radiative flux over the distance it has traveled inside 
the droplet, z/. 1 is found from Snell’s Law (neglecting the slight 


dependence of n on wavelength): 


n sin@' = sinO; uu =cos0'’; n=1.33_ 


Ht = cos [sin-! (0.75 sin@)] (2.2) 


E, 9 is found using I) ,, the spectral intensity at polar angle @ due to the 


radiative source, and the differential solid angle da: 


Exe =h,do where dw=dA, / R2 (2.3) 


dA,, only covers a fraction of the hemisphere at 6, that fraction of the 


hemisphere occupied by the source of the radiation. The author designates 
this fraction the fractional coverage with symbol fg, the f chosen because of the 
relation to the view factor nomenclature. fg will be found for the specific 


laboratory geometry using analytic geometry. 
Then dq = 2 n fgsin@ d@, and using the identity E = n I [16], we have 
F, 9(Z) = 2 Ey; fg cos® sin® (1 - pg) exp(-K,.z/p) dO (2.4) 


Integrating over all possible wavelengths and polar angles gives 


F(z) = 2 Ja Eas [saa fg cos® sin® (1 - pg) exp(-K,z/p) dO da (2.5) 


The volumetric heat generation due to absorption by radiation H appears 


because of conservation of the energy missing from F as z increases [17]: 


H(z) = -dF/dz (2.6) 
oo yi 
=2 IN Ex 5 Ky 1 (1/1) fg cos® sin® (1 - pg) exp(-K,z/p) dO da 


These important expressions for F and H are developed by Viskanta and Toor 
[17] for the more general case including scattering and bottom reflection. 

Their particular application was local absorption of solar radiation in a pond. 
Their expressions do reduce to (2.5) and (2.6) under the assumptions made in 


this development. 


In order to numerically evaluate the double integrals in the expressions for F 


and H, the following functions of 6 and A must be determined: 
(1) E,,9, the incident polar and spectral distribution 
(2) fg, the fractional coverage of the hemisphere as a function of 0 
(3) pg, the droplet surface reflectivity, a function of 0 
(4) xy, the droplet spectral absorption coefficient, a function of A 


The next four sections are devoted to this task. 


2.2 Assumption of Blackbody Source 
The three radiative heater coils of the experiment may be considered to be 
blackbodies. Then, the spectral distribution of the incident radiation is given 


by the familiar Planck distribution [16]: 


E,.6=Enp = Cy / {05 [exp (C2 / ATR) - 11) (2.7) 
where Cy =2 mh co? = 3.742E8 W-um4/m2 
Cy = (h cg / k) = 1.439E4 pm-K 


Of course, the radiative heater coils are not perfect blackbodies. However, 
this approximation is the only reasonable means of quantitatively 
characterizing the spectral dependence of the incident radiation. 
Furthermore, the integration step will likely smooth out any differences 


between reality and the blackbody assumption. The heater coil temperature 


TR is read as an experimental digital display. 


2.3 Analytic Calculation of Fractional Coverage fg 
With a few reasonable assumptions analytic geometry can yield the functional 
form of fg. These assumptions are listed in the following: 


(1) The axes of symmetry of the two upper cone-shaped heater coils pass 


(2) 


(3) 


(4) 


through the center of the liquid-solid interface. 

The elevations and radial positions of the two cone-shaped coils are 
equal. | 

The ring-shaped coil is horizontal and its axis of symmetry passes 
through the center of the liquid-solid interface. 


The small eccentricity of a given point on the droplet is neglected. 


Figure 3 is the vertical section of the geometry perpendicular to the faces of 


the cone-shaped coils with the needed measured dimensions. The 


measurements also implicitly assume that the two cone-shaped coils are 


radially opposite each other. The dimensions of Figure 4, a fabrication 


drawing of the cone-shaped coil assembly, are also required. 


Planar trigonometry gives the fg function corresponding to the ring-shaped 


geometry in a straightforward manner as illustrated in Figure 3. Note that fg 


is constrained to be Zero in a range of 6 because the lip of the ring-shaped coil 


blocks a section of each cone-shaped coil. The analysis of the two cone- 


shaped coils is more involved: 


(1) 
(2) 


(3) 


(4) 
(5) 


Find the center point of the front face of one of the cone-shaped coils. 
Use the radius of this disk-shaped face to find the equation of the 
sphere containing the circular edge of the disk. 

The intersection of this sphere and the sphere centered at the origin 
with radius equal to the distance from the origin to the point of highest 
elevation of the cone is the equation of the disk face of the cone-shaped 
coil. 

Write this equation of the circular disk face in spherical coordinates. 
Then calculate the azimuthal position of the disk $q as a function of 


polar angle 8. Because there are four total semicircular disk sections, 


cone-shaped 


heater \ 


ring-shaped 
heater 


center of wetted area 61 = 90° - arctan(51/131.5) = 68.8° 


6, = 90° - arctan(36/131.5) = 74.7° 
63 = 90° - arctan(10/131.5) = 85.7° 


263/2 = 131.5 mm 


FIGURE 3 
Measurements of Radiant Heater Panel Geometry 
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FIGURE 4 
Drawing of Conical Heater Panel 
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fo = 46q / 2n (2.8) 

(6) Subtract out the effect of the holes in the cone-shaped coils (truncated 
sections) by repeating the analysis using the center of the back face of 

the cone and the sphere centered at the origin and passing through the 


edge of the back disk face. 


The required nomenclature is presented in Figure 5 and the mathematical 


manipulations are given below. 


FIGURE 5 
Conical Heater Panel Coordinate System and Nomenclature 


(1) Choose o-.=0 
Then C = (f¢, 8c, 0) 
Xc = Pc SinOc cOSb¢ = Pc sinB> 
Yc = Pc sinO- sing, = 0 


Zc = Pc COSO¢ 
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(2) 


(3) 


(4) 


(5) 


(6) 


And (x - X¢)2 + (y - yo)? + (Z - 2)? = ra? 


Thus the equation for the sphere centered at C is 


x2 - 2xx¢e+ y2 + 22 - 2ZZ% = Ig? - Xe? - Ze? (2.9) 
The equation for the sphere centered at the origin is 

x2 + y2 +22 = p2 (2.10) 
Then 2xX¢ + 2ZZe = p2 - rg? + Xc2 + Zc? = 2p? 
Or xX == (Ze / Xo) Z+ Pc? / Xe (2.11) 


In spherical coordinates 


p sin® cosog = - (Ze / Xc) p cosO + Pc? / Xe 


Then ba = cos"! [p. / (p sin@¢ sin®) - cot, cot®] (2.12) 

And fgfront = (2/1) cos-1[p¢ / (p sin®@¢ sin®) - cot@, cot®] (2.13) 
where p = (pc2 + rg2)1/2 

Finally fg = fofront - fgback (2.14) 

where fgback = (2 / m) cos-1 [p,back/(pback sin@, sin) - cot®, coté] 


Of course fefront and fgback must be set to zero for certain ranges of 0. 


Using trigonometric identities one can show that the fg function does in fact 


become zero for the uppermost and lowermost elevations of the cone-shaped 


coil disk face. The arccosine argument becomes equal to one: 


First 04 BR =6.+7; Pc / Pp = cosy 
Then Pc / (p sin®- sin) - cot6, coté 
= [cosy / sin®, - cotO, (cos6, cosy F sin@, siny)] / sin(8¢. + ¥) 
= [cosy (1 - cos20,) + cos@, sin@, siny] / [sin®, sin(@, + y)] 
= (cosy sin26, + cos@¢ sin@¢siny) / (cosy sin2@, + cos@¢ sin@,-siny) 


=1 
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Figure 6 shows the functional forms for the cones and holes. Note that the 
functional form is biased toward smaller polar angles; less surface area need 
be covered at higher elevations. Figure 7 shows the complete fg function after 
the subtraction of the holes and addition of the ring-shaped coil. The function 
was saved numerically in a data file for 0.1 degree increments for use in the 
numerical double integrations of F(z) and H(z). Note that fg is fairly sensitive 


to the geometrical measurements Pp, and 6.. 


2.4 Reflectivity of Liquid Surface pg 

Because water can be considered a dielectric, its reflectivity is strongly a 
function of the incident angle (already assumed to equal the polar angle 98). 
By assuming that the incident radiation is randomly polarized the reflectivity 
Pg is found as the arithmetic average of the perpendicular and parallel 
interface reflectivities given by the Fresnel relations of electromagnetic theory 


[18): 
Pe = (Ri + Ry) / 2 (2.15) 
Because water is weakly absorbing (k << n), the effect of the absorption 


optical constant k, and its dependence on wavelength, can be neglected. Then 


the Fresnel relations reduce to 


R, = [(cos@ - u) / (cos@ + u)]}2 (2.16) 
Ry = [(n2 cos® - u) / (n2 cos@ + u)]? (2.17) 
where u = (n2 - sin2@)1/2 


and again n = 1.33 is a sufficiently accurate value representative of all 
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FIGURE 6 

Fractional Area Coverage of Cone Face and Hole 
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FIGURE 7 
Fractional Area Coverage fg 
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wavelengths. The plot of pg in Figure 8 shows that radiation incident at 
shallow angles will not enter the droplet-solid system but rather will be 


reflected back to the surroundings. 


2.5 Absorption Coefficient Kk, 

The absorption coefficient x, is just the exponential damping coefficient of 
the radiative energy flux with distance through the absorbing medium. For 
water k, is a very strong function of wavelength with water generally 
appearing opaque to large-wavelength, low-energy radiation and transparent 
to small-wavelength, high energy radiation. Because the blackbody 
distribution is also a strong function of wavelength, it is prudent to use a 
relatively precise form of the x, function from the literature. Reference [19] 
gives the optical constant k for many wavelength intervals, and that data is 
used for the model presented here. The electromagnetic theory provides the 


conversion to «} [18): 


Ky, =4nk/A (2.18) 


Care was taken to accurately characterize the function for large wavelengths 
(> 10 tm), because a significant percentage of the energy is found in this range 
for the radiative heater temperatures available in the laboratory setup. 

Figure 9 is a log-log plot of k, versus wavelength. Planck blackbody 
distributions are superimposed on the plot and approximate the spectral 


dependence of the incident energy for several values of Tr. 


2.6 Numerical Integration Scheme 


The numerical double integrations for F(z) and H(z) were carried out using a 


i ee) 


Reflectivity rho 


Azimuthal Angle thetaO 


FIGURE 8 
Reflectivity pe vs. Polar Angle @ 
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Spectral Absorption Coefficient of Water and Plank Blackbody Distributions 
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simple rectangular integration scheme, but the fineness of the integration grid 
guarantees sufficiently precise results. The step size in 6 was 0.1 degree while 
the step size in 4 varied to correspond with the optical constant data. The 
analytic expressions for E) p, 1, and pg placed no additional constraints on 
step size. Integration over wavelength was carried out to 200 um, where the 
incident radiation can certainly be neglected (E, , << 0.1 W/m2). The lower 
limit on wavelength was chosen to be 0.2 um where again the incident 
radiation can be neglected, but here only for radiative source temperatures 
up to about 3000 K (the small bandwidth from 0.0 to 0.2 m justifies 
neglecting this energy for higher temperatures). Appendix A presents a 
listing of the QuickBASIC computer code used in the actual computations of 


F(z) and H(z). 


2.7 Results 

Figures 10 through 12 graphically display the calculation results for the F(z) 
function, the radiative flux at a depth z below the liquid-vapor interface. The 
normalized quantities are with respect to the value at the liquid-vapor 
interface using the corresponding radiative heater temperature. Figures 13 
through 16 display the results for the H(z) function, the volumetric heat 
generation term, using the same normalization scheme. Figure 16 is included 
only to show that H(z) truly becomes negligible for low values of Tr. The 


laboratory range of radiative source temperatures Tp is from 475 to 650 C (750 


to £25 K) while 5800 K is the characteristic solar temperature. 
Several interesting and important conclusions can be drawn from the graphs: 
(1) For laboratory Tp, the radiative flux drops off to approximately 7% of 


its surface value at a depth of 0.1 mm. 
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FIGURE 10 
Radiative Heat Flux F(z) (z < 2 mm) 
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Normalized Frad(z, Tcoil)/Fradmax(Tcoil) 


Depth z (mm) 
FIGURE 11 
Normalized Radiative Heat Flux F(z) (z < 5 mm) 
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FIGURE 12 
Normalized Radiative Heat Flux F(z) (z < 0.2 mm) 
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FIGURE 13 
Volumetric Radiative Heat Generation H(z) (z < 0.2 mm) 


veccececsccectsceccecscoscscbcscencecccccsdececosencncesgecsscecsscsaceherstecereasee docenceasessreepeccnscenserseqenncnscsnssesepoeccscsaeans 


Normalized Hrad(z, Tcoil)/Hradmax(Tcoil) 


- : : : 

obec ccc cc ccc scence pe ccccs ccs coscqcescnnsvessesegs saenssencess qosccccccescscegncseesescorseusessessacsssanpesserscessuns 
cere . 3 . * * * ’ 
a ‘ * ‘ ' 


0.00 — 
0.000 0.020 0.040 0.060 0.080 0.100 0.120 0.140 0.160 0.180 0.200 


FIGURE 14 cre) 
Normalized Volumetric Heat Generation H(z) (z < 0.2 mm) 
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Normalized Hrad(z, Tcoil)/Hradmax 
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FIGURE 15 
Normalized Volumetric Heat Generation H(z) (z < 0.2 mm) 
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FIGURE 16 Peraaen mee) 
Normalized Volumetric Heat Generation H(z) (z < 2 mm) 
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(2) | The volumetric generation term drops off to approximately 6% at a 
depth of 0.02 mm. 

(3) For laboratory Tp, the volumetric generation term drops off to 
approximately 0.4% at a depth of only 0.1 mm. 

(4) For higher Trp, the volumetric generation term becomes roughly linear 
(and to a lesser approximation constant) for z greater than about 


0.2 mm. 


These observations suggest that the radiative energy reaching the droplet may 
be input to the droplet-solid system at three locations: the liquid-vapor 
interface, the droplet interior, and the solid-liquid interface. In addition, for 
laboratory conditions, the term input to the droplet interior may be 
neglected. The details of the simplified model are given in the following 


subsection. 


2.8 Simplified Model 

Because the method of adding the volumetric heat generation term must be 
compatible with the liquid heat equation models used in the overall 
computer model, it is unrealistic to use the entire functional form of H(z, Tp). 
The typical droplet has an initial apex of from one to four millimeters, so any 
reasonable finite difference grid could not represent the H(z) curve. The 
equation for the analytic quasi-steady state liquid solution, given in Chapter 
3, requires two integrations of the H(z) function. These analytic integrations 
are not practical because of the difficult curve fit required (the next 
subsection only gives a piecewise curve fit). The solution is to take advantage 
of the very thin layer in which much of the absorption takes place and assume 


a surface boundary condition and perhaps a residual H term. 
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The radiative flux that penetrates into the droplet has a vertical value just 
inside the droplet given by (2.5) with z = 0. This flux is calculated with help 
of the quantity F, as defined by the author: 


oo 2 
F(z=0") =? J, Ey b like fg cos8 sin®@ (1 = Po) dé da 
oo Z 
= U Exb dA] i if fy cos@ siné (1 - Po) de] 


= Ep F, 
= F, oO Tr? (2.19) 


where numerical integration gives 
¢ = 0.2261 (2.20) 


and the Stefan-Boltzmann law has been used. 


Because H(z) drops off so quickly just inside the droplet, the flux given by 
(2.19) can be divided into three portions: that absorbed at the Rapiahoapien 
interface, that absorbed throughout the droplet depth, and that absorbed at 
the solid-liquid interface. Assuming as before no reflection at the droplet 
bottom, the portion absorbed at the solid-liquid interface is a fraction of the 


quantity given in (2.19) to be found by setting z to 6 in (2.5): 


F(z=8) = F, 6 Tp4 (FRACTION at z = 8) (2.21) 


This term is subtracted from the liquid-solution conductive flux for use by 
the solid BEM and will become significant when the droplet becomes 
extremely thin near the end of the evaporative process. The next subsection 
gives a curve fit for (FRACTION at z = 8) so that the code can use (2.21) 


without storing a large numerical array of the radiative flux as a function of z 


Pata 


and Tr. 


Figure 15 suggests that a uniform value for the residual volumetric heat 
generation H, is a reasonable approximation. Then, the liquid-vapor 
boundary condition flux is calculated using an energy balance (per unit 


horizontal area): 


F(z=0) = F, 0 Tp? - F(z=8) - He 8 
= F, o Tr4 [1 - (FRACTION at z = 8)] - He 8 (2.22) 


For energy input to the droplet interior, the He term is also used in a 
nonhomogeneous term in the tridiagonal finite difference solution and in the 
analytic quasi-steady state solution. Note that H, is applied in the entire 


droplet thickness including the thin primary absorption zone at the liquid- 


vapor interface. A simple curve fit for an H<(Tp) function could be carried 
out. Because H,; is extremely small for radiative source temperatures under 
laboratory conditions, the EVAP code itself assumes that H, is zero. 

However, the development of Chapter 3 only assumes that He is uniform, so 
the EVAP code may easily be extended to cover higher values of Tr. Rather 
than assuming a constant form for the residual heat generation, a linear profile 
could be used. However, this would significantly complicate the quasi- 


steady state analytic solution, and this option is not explored. 


2.9 Curve Fit for Radiative Flux F(z, Tp) for Laboratory Tp 
The curve fit developed here is used in the overall model FORTRAN code in 
order to quantify the radiation directly absorbed by the solid in the wetted 


region. Figure 17 suggests a hyperbolic or exponential form for a curve fit. 
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FIGURE 17 
Normalized Radiative Heat Flux F(z) for Laboratory TR 
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FIGURE 18 
Inverse Normalized Radiative Heat Flux F(z) for Laboratory Tr 
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The author chooses to try a hyperbolic-type fit by plotting the reciprocal of 
normalized F(z) shifted down one unit (Figure 18). The six curves, 
representing the laboratory Tp temperatures, are indeed somewhat linear. 
Piecewise linear fits are made over the 1/F(z) - 1 curves over five ranges for z 
less than 1 mm using least-squares statistical spreadsheet functions. The range 
for Z up to 1 mm is also used for z greater than 1 mm because the curves 
become rather linear at z = 1 mm and because normalized F(z) becomes rather 
small for z greater than 1 mm. Each linear fit for a range and Tp value results 
in a slope m and y-intercept b. Another spreadsheet program allowed easy 


least-squares second and third order polynomial fits of the m and b 


quantities. The resulting form of the curve fit of F(z, Tp) is 


F(z, Tr) = F, 6 Tr4 / (m(Tp)z + b(TR) + 1) (2.23) 
(FRACTION at z = 6) = 1 / (m(TR)d + b(TR) + 1) (2.24) 
where m(TR) and b(Tp) are listed in the final FORTRAN code of 
Appendix B in the FUNCTION FRAD(D, T). 


The curve fit for H(z, Tp) is easily found by partial differentiation: 
H(z, Tp) = OF/dz = F, 6 Tr4 m(TR) / (m(TR)z + b(TR) + 1)2 (2.25) 


Figure 19 shows that the original F(z) function and its curve fit evaluated for 


TR of 500°C closely match with the maximum error being 3% of the liquid- 


vapor interfacial value. Figure 20 shows that using the curve fit in (2.25) to 


calculate H(z) degrades the agreement somewhat. 
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Frad(z, Tcoil=500C) kW/nm2 


FIGURE 19 Sea) 
Comparison of F(z) and Curve Fit for Tr = 500°C 
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FIGURE 20 Beri een 


Comparison of H(z) and Curve Fit for TR = 500°C 
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2.10 Summary of Assumptions and Errors 

The assumptions of this chapter, both in the governing equation (2.6) and in 

the simplified model, have made the problem of direct radiation absorption 

by the droplet tenable. The assumptions and corresponding justifications 

used to arrive at the radiation model used in the computer calculations of 

Chapter 4 are summarized in the following: 

(1) blackbody behavior of the radiation emitting heat panels—panels 
have black coatings 

(2) geometrical symmetry including azimuthal symmetry of the droplet 
and symmetry of the radiative panels with respect to the droplet— 
droplet symmetry before impact and experimental design 

(3) reflection at the liquid-vapor interface obeys the Fresnel relations for 
randomly polarized incident radiation—water can be considered a 
dielectric 

(4) the liquid-vapor interface is horizontal and flat—droplet rather thin 

(5) scattering of radiation within the droplet is negligible—droplet is 
relatively thin, and the water is degassed 

(6) no reflection at the solid-liquid interface—emissivity of solid 
approximately 0.84 

(7) radiation entering the droplet is completely absorbed at the surface— 


justified by evaluation of the very small residual volumetric heat 


generation for laboratory values of Tp. 


Of the seven the most questionable assumption is (4) because the outward 
normal near the edge of the drop may in fact be closer to horizontal than 
vertical. Equation (2.6) then overestimates the fraction of energy reflected 


away from the droplet. The conclusion that under laboratory conditions 
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absorption is restricted to an extremely thin zone at the surface would not be 

affected, but two options for improvement in the overall code (Chapter 3) do 

become apparent: 

(a) | Assume some value for F, between that calculated with assumption (4) 
(i.e., F, = 0.226) and that resulting with pg = 0 (i.e., F, = 0.255). 

(b) Make F, a function of the local geometry through r and the derivative 
(slope) of the droplet profile. 

The relatively small 13% difference between the two F, values suggests that 

using assumption (4) is not critical to the simplified radiation model so long 

a: Ig may be neglected although this hypothesis has been put forward to 


explain the results of Chapter 4. 


On the other hand, for real-fire non-laboratory conditions the removal of 
assumption (4) may affect the value of H, and make it a function of r to be 
evaluated as the droplet profile is updated. These evaluations could be 
expected to require significant and intensive computational effort. Not only 
would the local reflectivity affect H, but so would the variable distances to 


the free surface (lengths for exponential decay). 
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3. SECTION II - COMPUTER CODE DEVELOPMENT 

This chapter reviews in detail the development of the overall model 
FORTRAN code (PROGRAM EVAP) listed in Appendix B. Although the 
code was completely re-written, previous workers in this NIST-Center-for- 
Fire-Research-sponsored line of research have written codes with many 
similar elements [8, 10-12], and EVAP builds on their work. Re-writing 
allowed for maximum understanding of the effect of each level of model 
refinement while comparison with previous models and results provided a 


firm base of reference. The following program features are new: 


(1) radiation model of Chapter 2 incorporating the experimental geometry 
of Dawson [13] 

(2) framework for radiation model allowing for constant residual 
volumetric heat generation term—needed for high Tp 

(3) exclusively one-dimensional liquid model 

(4) initial liquid temperature distribution given by the error function 
solution to the semi-infinite solid problem 

(5) | avoidance of computationally intensive, extremely small initial time 
steps by assignment of initial quasi-steady state status to outer droplet 
nodes based on Fourier number criterion 

(6) smooth transition from transient tridiagonal finite-difference solution 
to quasi-steady state analytic solution 

(7) use of Clausius-Clapeyron equation in analytic linearization of 
evaporative liquid-vapor boundary condition 

(8) option of two droplet geometry models with the second model 


(subroutine written by S. Tinker) allowing detailed consideration of the 
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transient droplet contact angle @ and the receding angle 0, 

(9) detailed reporting of meaningful transient quantities such as area- 
averaged upper and lower droplet surface temperatures 

(10) criterion for termination of boundary element method (BEM) 


integration in recollection time 


Furthermore, the main code, subroutines, and functions are well annotated 
with all variables explicitly defined, and those quantities for which 
modification requires additional changes in the code are clearly identified. 
Effort has been made to make the 1500 lines of code readable and 
computationally efficient. Early versions of the program were written in 
QuickBASIC and run on PCs with 486 microprocessors. This allowed for 
sophisticated debugging capabilities and computational speed. Later 
versions were edited using the 486s and transferred to the NIST Tiber system 
for execution in FORTRAN so that the required bessel and error functions 
could be called and the code run with acceptable speed. The final version of 


the code requires approximately 30 minutes to run on Tiber. 


Because of high solid temperature gradients, a Green's function solution 
equation or BEM method is used to calculate the solid surface temperatures. 
The next subsection on the BEM precedes the subsection on program 
structure because of the strong influence BEM exerts on the structure. 
Subsection 3.3 describes the solid-vapor boundary conditions while 
subsections 3.4 through 3.9 are needed for the much more complicated solid- 
liquid boundary condition, i.e. the liquid solution. Subsection 3.10 shows 
how the BEM and liquid solutions are legitimately coupled. Finally, the last 


subsection details how the many property values were specified. 
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3.1 The Boundary Element Method (BEM) 

Because EVAP is required to model solids with low thermal conductivity, the 
spatial and temporal temperature distributions in the solid are crucial to the 
evaporative process. DiMarzo [10] noted results that finite-difference 
methods within the solid with reasonable time steps fail because of sharp 
localized gradients at the droplet edge. However, a boundary element 
method is well suited to the problem. By linearizing the solid-vapor 
boundary condition, the solid problem (governing heat equation and 
boundary conditions) becomes linear and a Green's function approach is 


valid. 


A heat conduction Green's function is a solution to the heat conduction 
problem having the same geometry but having homogeneous boundary 
conditions as the original heat conduction problem. It represents the 
temperature response at a point in time due to a heat pulse at another point at 
some previous time. By integrating all the heat pulses over space and 
recollection time, the current temperature distribution is found. Beck, et al., 
[20] present a very good full-text treatment of heat conduction using Green's 
functions in which they derive the general Green's function solution equation 
(GFSE) for heat conduction using mathematical manipulations such as 


Green's theorem. The general GFSE is 
T(r, t) = Tinit(t, t) + Tgen(r, t) rt: Tbe(r, t) (3.1) 
where Tinit is the temperature contribution due to the initial conditions, Tgen 


is the contribution due to volumetric heat generation within the solid, and 


Thc is the contribution due to the boundary conditions. Tgen is Zero because 


SP 


their is no internal heat generation, and Tinit can be conveniently be made zero 


by subtracting away the initial temperature distribution. 


The mathematical statement of the problem in the semi-infinite solid domain 


(see Figure 2) is 


dT/dt = as V2T (3.2) 
t = 0: T =Tsg-qoz/k, (3.3) 
O<r<R;z=0: -k, OT/dz given (3.4) 
r>R;z=0: -ks OT/0Z = ho (Ts - Ta) - Feo Tr? (3.5) 
Z— —oo; all r: -ks OT /0z = qo (3.6) 
r— 00; Z <0: OdT/or = 0 or equivalently -ks dT /dz = qo (3.7) 


Note that qg is negative for the case of radiative heat input. Tinit is made zero 


by mapping the problem to the u domain: 


u(r, Z, t) = T(r, z, t)- Ts9 + qQgz/k, (3.8) 
du/dt = a, V2u (3.9) 
t = 0: u=0 (3.10) 
O<srsRz=0: -k, du/dz given (second kind) (3.11) 


r>R;z=0: -ks du/dz = hg (Ts - Ta) - FE 6 Tr? - qo (third kind) (3.12) 
Z — —»o; all r: du/dz = 0 (3.13) 
r— co; 7 <0: du/dz = 0 (3.14) 


Only the boundary condition on the solid surface remains nonhomogeneous. 
Now using the GFSE with the expression for Tpc (second and third kind) 


yields an expression for the solid temperature distribution: 


aS 


T(r, z, t) = u(r, z, t) + Tsg- qQgZ / kg (3.15) 
- ou/odz G(r, z, t; r', 0, t') 2xr' dr' dt’ 


t ) 
where u(r,z,t)=a5J J 
t=O0r = 


and the proper form of the Green's function (typically found using Laplace 


transforms, separation of variables, etc.) is [20] 


GO szati Hd: t) 
= 2 (4nagt')"3/2 exp(-z2/4ast') exp[-(r-r')?/4ast'] Lo(rr'/2act') (3.16) 


where Lo(a) = e-lg(a) is the exponentially damped modified Bessel function of 
the first kind of order zero and is calculated by calling a system FORTRAN 


subroutine. 


Note that t', the recollection time, is zero at the current time step and increases 
in the negative time direction, i.e. t' = t - tg where tg is a dummy integration 
variable in the forward time direction. The liquid solution and infrared 
thermography experimental validation only require the solid surface 
temperatures, so the simplification z = 0 is made.in the resulting expression 


for u: 


u(r, z=0, t) (3.17) 
= (4na,)"1/ 2h sate p(ou/az) r’ t'3/2 Lo(rr'/2ast') exp[-(r-r')2/4at'] dr’ dt' 


where du/dz = 8T/dz + Gg / k, = (qg- qs) / ky (3.18) 


The numerical evaluation of (3.17) at each time step is referred to as the BEM. 


Two difficulties become apparent: a) at what r do the integration terms 
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become negligible? and b) how to handle the degenerate cases when r or r' 
may be zero or equal or when t' is zero. Liao [12] showed that r = 10Rg is 
sufficiently large for an upper bound. Zero points are avoided by a 
nodalization scheme in which functions are evaluated in the center of radial 
and time steps. The case for integration terms where r = r' is handled by an 
approximate analytic integration just over the element of concern in which an 
exp(u2)-type term becomes an erf-type term. The resulting numerical 


integration expression as derived by Liao [12] and verified by the author is 


ur, Z=0,.t) (3.19) 
= (4n0,)"1/ De Z_  (u/dz) r t'-3/2 Lo(rr'/2ast') exp[-(r-r')2/4arst'] Ar’ At' 


os y, (du/dz) r t'-! Lo(r2/2ast') erf [0.5Ar/(4act')1/2] At’ 


where T(r, z=0, t) = To(r, t) = Tog + u(r, z=0, t) (3.20) 


The temporal integration step size must be chosen sufficiently small to 
capture the temporal behaviors of the Green's function and u gradient while 
the overall, forward-marching time step should only capture the temporal 
behaviors of the surface fluxes and temperatures. The author chooses overall 
time steps of 0.1 s and 1.0 s with the smaller time step to be used for several 
seconds after droplet impact when the heat flux rapidly decreases according 
to the semi-infinite solution. These values are consistent with the time step 
chosen by Liao [12]. The products of the terms in (3.20) excluding the du/dz 
forcing function are considered the weight functions. Figure 21 shows the t' 
behavior of several typical limiting case weight functions (limiting values of r 
and r') where the weights which increase as t' approaches 0 apply for the case 
r=r'. The 0.1 and 1.0s time steps will not adequately capture the Green's 


function. EVAP handles this by integrating the Green's function itself over 
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FIGURE 21 
Typical Limiting-Case Weight Functions vs. Recollection Time 


the first second or two using a step size of 0.004 s and saving the values in a 
weight matrix. For larger recollection times u is calculated at the center of the 


overall (0.1 and 1.0 s) time steps. 


The change from the 0.1 to 1.0 s time step requires careful handling of the 
matrix of past forcing functions, the past solid surface heat fluxes mapped 
into the u domain by subtracting away the initial uniform flux. When t equals 
4s, the time step At is changed to 1.0 s, and SUBROUTINE RECONF is called. 
The algorithm in RECONF redefines the forcing function matrix FRCFNC so 
that the forcing function values at the centers of the 1.0s time steps are used. 


The key line of code is 
FRCFNC(K, J) = (FRCFNC(10*K - 5), J) + FRCFNC(10*K - 4), J)) / 2 


Because the Green's function decays with recollection time, integration all the 
way back to droplet deposition may not be necessary. However, when the 
droplet becomes completely evaporated the required maximum recollection 
time or "memory" changes drastically. EVAP uses a dynamic measure of the 
required "memory" by aborting the time step integration loop if all the ratios 
of the current u addition to the maximum u addition (all nodes) are less than 
a user-defined ratio called NEGRAT. Results for the output MEMORY are 


given in Chapter 4. 


The radial step size becomes larger for large r where temperatures and heat 
fluxes are known to be much more uniform. The number of elements is rather 
conservative at 78 with 12 making up the initial wetted area, and the 


outermost element extends to a distance of 13Rg. The breakdown of element 
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size is given in FUNCTION WGHT. 


The steps that EVAP uses to calculate the current solid surface temperature 

given the vector of current solid surface fluxes (u gradients) are 

(1) Calculate and save integrated Green's function values for small 
recollection times in matrix called W. 

(2) | Update the matrix of past forcing functions FRCFNC using current 
solid surface fluxes. 

(3) | Call SUBROUTINE BEM1 (or BEM2) to calculate u vector as a matrix 
multiplication. 

(4) Call FUNCTION WGHT to calculate the needed weight (effect of one 
element on another or itself for a certain recollection interval). 

(5) For short recollection times (<1.0 s for At = 0.1 s and <2.0 s for At = 1.0 s) 
get the weight from W, and for longer recollection times explicitly 
calculate the weight using one formula for r #r' and another for r =r’. 

(6) Following completion of the multiplication of the weight and forcing 
function matrices, add u to Tgp to get Ts. 

These functions of EVAP were validated by assuming a constant and uniform 

flux over a circular disk on the semi-infinite solid, a problem for which an 

exact solution is known. The results of this validation are presented in 


subsection 4.8. 


3.2 Program Structure 

With the requirements of the solid solution (BEM) part of the code clear, the 

required structure of the remainder of the code is evident. The remainder of 
the code must take the solid surface temperatures from the BEM and return 


the solid surface heat fluxes. For the non-wetted region this is a simple task, 
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but for the wetted area a liquid solution is needed. 


The liquid solution assumes that all heat flux is one-dimensional with no 
radial component (previously shown to be small [15]), and proceeds through 
three modes: a) error function solution based on the instantaneous contact of 
two semi-infinite bodies; b) tridiagonal transient solution; c) quasi-steady 
State analytic solution. Boundary conditions consider the effects of radiative 
input, convection, evaporation, and re-radiation at the liquid-vapor, solid- 
vapor, and solid-liquid interfaces. The evaporative boundary condition 
requires careful calculation using the Chilton-Colburn analogy and Clausius- 
Clapeyron equation while the radiative boundary condition relies on the 
simplified model from Chapter 2. Finally, because of the critical and 
unpredictable effect of droplet impact, surface tension, and flow, two 
empirical models are used to specify the transient droplet profile. Finally, 
the overall model must couple the liquid and solid domains using a simple 
predict-correct method in which there is feedback between solid surface 
temperature and heat flux. Figure 22 is a simplified flow chart of the program 


structure. 


3.3 Solid-Vapor Boundary Condition 

The non-wetted solid surface boundary condition must account for radiation 
absorbed from the heat panels, re-radiation from the solid surface to the 
surroundings, and convective cooling to the ambient. The boundary 


condition on T(r, z, t) can then be written: 


T> R; Z= 0: -ks oT /dz = Gconv + dre-rad = Grad (3.21) 


So 
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FIGURE 22 
Simplified Code Flowchart 
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FIGURE 22 (Continued) 
Simplified Code Flowchart 


The convective term can be written using the experimentally measured 
convective heat transfer coefficient [21] as a function of initial solid surface 


temperature Toco, curve fit by the author (see Figure 23): 


Qconv = h (Ts - Ta) where 
h = -42.348 + 1.3663Tcg - 0.011498Tc92 + 3.1954E-05T 59? W/m2-K (3.22) 


with 70°C < Tso < 160°C 


The linearized form of the re-radiation term is given because the BEM 
requires linear boundary conditions and because the solid surface 


temperature is expected not to vary more than about 50°C: 


Gre-rad = h, (Ts - Ta) where 
he =e ((len- 15 Chee Glaneabneiter 1.2) PWGER 


and (Tso - 15°C) is a reasonable value of the expected solid surface 


temperature in order to make h, constant. The emissivity of the Macor tile is 


given as € = 0.84. The radiation absorption term is theoretically calculated 


using the fg term from Chapter 2: 


oo /2 
Grad = 2€ J. Exp! ba fg cos® sin8 d@ da 
oo 2 
=e U. Ep dA) Ue 2 fg cos® sin® dé] 


=€E Ep F 
=FeoTp4 (3.24) 
where numerical integration gives § F = 0.2552 (3.25) 


By defining the overall heat transfer coefficient as ho = h + hy, one arrives at 
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Experimental Convective Heat Transfer Coefficient vs. Tsg (Data from [21]) 
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equation (3.5), the solid-vapor heat flux needed by those boundary elements 


not covered by liquid: 


r>R;z=0: -k, OT/0z = hy (Tg - Ta) - Feo Tr4 (3.5) 


3.4 Semi-Infinite Liquid Solution 

The purpose of the semi-infinite liquid solution is to handle the singularity of 
the infinite initial heat flux caused by the temperature discontinuity between 
the solid and just-deposited droplet. The time step required by the 
tridiagonal solution of the next subsection to accurately handle the high heat 
fluxes would be too small to be practical. However, before the diffusive heat 
wave reaches the liquid-vapor interface, the liquid temperature distribution 
will (to an acceptable approximation) be that of a semi-infinite body brought 
into instantaneous contact with another semi-infinite body, the solid. Another 
needed but reasonable assumption is that the interfacial solid-liquid solution 
is constant for the short time that the heat wave has not yet reached the liquid- 
vapor interface. Then the classic solution by Carslaw and Jaeger for a semi- 
infinite body initially at a uniform temperature with constant surface 


temperature can be used: 


T(r, z, t) = Te - (Te - Ty) erf (0.52z(ay t)-1/2) (3.26) 
qs = ki (Tc - T1) (n ay t)-1/2 where (3.27) 
Te = [(kipicy)1/2T + (kspscs)!/2Tso] / [(kipicy)1/2 + (kspscs)!/2] (3.28) 


is derived by equating the two semi-infinite body surface fluxes. Note that 
(3.28) is strictly valid only for a constant not linear initial solid temperature. 


However, the heat flux from the initial linear temperature distribution is only 


a small fraction of the heat flux given by (3.27), and neglecting this effect is 


certainly within the level of approximations already made. 


This solution is enforced until the heat wave reaches the liquid-vapor 
interface, arbitrarily assigned to be when the error function becomes about 
99.6% of its maximum value of 1.0. This percentage was chosen because it 
means that a 1.0 mm thick region of the droplet will require 0.4 s before the 
heat wave hits the upper surface. Then the time for switching to the 


tridiagonal solution tsjs is calculated as 


erf (0.55(c4 tsts)-1/2) = 0.9962 
0.58(cy tsts)-1/2 = 2.0611 
tsis = 62/ 17Q) (3.29) 


In addition, the radiative, convective, and evaporative boundary conditions 
need not and cannot be applied at the liquid-vapor interface during the semi- 
infinite solution because conditions there have not changed since before 
deposition (t = 0) and are considered to be off at infinity. The time tsIs varies 
for each wetted solid surface element because the thickness of water above 
each element is different. The initial liquid temperature profile for the 
tridiagonal solution is found using (3.26) except that Tc is replaced by the 
current temperature calculated by the BEM in order to keep the profile 


smooth as required by the tridiagonal solution for physically meaningful 


results. The BEM cannot function with the restriction that Ts be constant. 


3.5 Tridiagonal Transient Liquid Solution 


The droplet liquid can be assumed to be stagnant, and the small Eulerian 
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temperature time derivative associated with droplet edge recession 
(geometry Model B) can be neglected. Then conduction and radiation 
absorption are the only modes thermal energy transfer. Furthermore, the 
radial conduction term is small and can be neglected [15]. The resulting one- 
dimensional transient heat equation with volumetric heat generation and 


applicable boundary conditions are 


OT /dt = a 02T/dz2 +a He /k (3.30) 
T(z=0) = To(t) (3.31) 
OT(z=5)/d0z = A T(z=8) + B (3.32) 


where the solid-liquid surface temperature is given by the BEM, and the 
liquid-vapor boundary condition has been linearized as is required to make 
the liquid solution linear and amenable to the tridiagonal method. 


Expressions for A and B will be derived in subsection 3.7. 


This transient heat problem is numerically solved by dividing the droplet 
into annular columns with each column atop a BEM element. Each column is 
divided into n-2 finite difference elements with a phantom node just across 
each boundary. This arrangement with no nodes on the boundaries better 
represents the nonhomogeneous volumetric heat generation terms. Then the 


discrete form of (3.30—3.32) is 


(Tj* - Ti) / At = @ [(Tj-1* - 2T)* + Ti41*) + (Ti-1 - 27) + Ti41)] / 2 (Az)2 


+aH./k for i = 2, 3,..., n-2, n-1 (3.33) 
(T}+ +.T>*)'/ 2= Ts (3.34) 
(Txt by Tn-1*) / AZz=A (Ty-1* + Tn*) / 2+B (3.35) 


where the second order central derivatives are averaged to take advantage of 
the guaranteed stability of the Crank-Nicholson method. The small 

temperature changes caused by the very small velocities of the nodes as they 
compress together as the droplet gets smaller can also be neglected. Moving 


the unknown new temperatures to the left in this implicit system yields 


-y Ti-1t + (1 + 2y) Tit - y Tis = Ti + y (Ti-1 - 2Tj + Ti+1) + 7] 


fori =2,3,...,n-2, n-1 (3.36) 
0.5T y+ + 0.5T2t = Ts | (3.37) 
(-1 - 0.5AAz) Tp-1* + (1 - 0.5AAz) Tht = B Az (3.38) 
where y=aAt/2(Az)2 and n=aH,At/k (3.39) 


An efficient standard routine, part of SUBROUTINE GAUSEL, is used by 
EVAP to solve this system by arranging the coefficients in an nx4 tridiagonal 
matrix. The residual heat generation H, need not be taken to be unifrom in 
(3.39). Because H¢- was found to be negligible in Chapter 2 for laboratory 


conditions, it is set to zero in the code. 


3.6 Steady-State Liquid Solution 

Although the Crank-Nicholson method is unconditionally stable, its 
numerical accuracy is dependent upon the value of the nondimensional 
multiplier y= «a At / 2 (Az)?. Al-Khafaaji and Tooley [22] recommend that y 
be limited to 0.5. Unfortunately, as the droplet shrinks the axial step size Az 
causes ¥ to quickly increase. However, during this time, the temperature 
distribution in each annular column becomes smoothed out, and the time 
derivative dT/dt becomes negligible (as will be documented in Chapter 4). 


Then the liquid problem can be solved analytically yielding expressions for 
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the solid-liquid heat flux and updated liquid-vapor interfacial temperature. 
The current liquid-vapor interfacial temperature is needed to evaluate the A 
and B quantities in the linearized boundary condition and the evaporative 
mass flux. This solution is referred to as quasi-steady state because the 
boundary conditions and droplet thickness can still vary with time but not 
fast enough to warrant the need for evaluating dT/dt. The quasi-steady state 


problem statement and solution follows: 


d2T/dz? = -H,/k (3.40) 
T(z=0) = Ts (3.41) 
dT;/dz = AT; +B where here T; = T(z=8) (3.42) 
T(z) = -H,z2/2k + c1z+ 2 and dT(z)/dz = -Hoz/k +c 

c2 = T(0) =Ts 

T; = -H,.62/2k + c18 + c2 and dT;/dz = -H,8/k + cy 

then, -H,8/k + cy = A(-H,52/2k + c75 + co) + B 


solving for ci gives 


1 = (AH(52/2k - Hp5/k - ATs - B) / (A8- 1) 


qs = k (AH(52/2k - H¢5/k - ATs - B) / (1 - A8) (3.43) 
T; = -H¢52/2k - (AH,82/2k - H¢d/k - ATs - B)8/(1- A8)+Ts (3.44) 


Note that H, must be constant in the above derivation. The effect of He is to 
bow out the radiation-free linear temperature solution thus pushing thermal 
energy toward both interfaces. More complex forms of He may be used, but 
the resulting derivation would be somewhat complicated. Again H, is small 


enough to be neglected for the laboratory conditions modelled by EVAP, so 
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the following expressions are substituted for (3.43) and (3.44): 


qs =-k (A Ts + B) / (1- A8) (3.45) 
T; = 5(A Ts +B) / (1- A 8) + Tg = (Tg + B8) / (1-A 8) (3.46) 


3.7 Liquid-Vapor Boundary Condition and Volume Flux 
Thermal energy is transferred at the liquid-vapor interface by radiation 
absorption (20 zm deep layer), convective cooling, and evaporation. 


Therefore the liquid-vapor boundary condition is 

-k 0T;/0Z = qcond = Aconv - Grad + Gevap (3.47) 
Radiation 
The radiative flux assumed to be input at the surface is given by (2.22) and 


(2.24): 


Grad = F, 6 Tr4 [1 - 1 / (m(TR)5 + b(TR) + 1)) - He 8 (3.48) 
Convection 
The droplet is assumed not to disturb the natural convection patterns 


established before deposition. In addition, the convective term will be seen to 


be only 1-2% of the total flux. Then the convective flux is 


Gconv = h (T; - Ta) (3.49) 


where h is given by (3.22). 
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Evaporation 
The Chilton-Colburn analogy and Lewis number Le predict a mass transfer 


coefficient h,, which when coupled with the expression for total molar flux N 


will lead to an expression for evap: 


N} = [diffusive flux] + [bulk motion] = J)* + C) v* (3.50) 
using Fick's Law, Nj =-CD dxi/dz + x} (Ni + Nair) (3.51) 
neglecting the solubility of air in the droplet, Nair = 0 

then N] = -CD (dx;/dz) / (1 - x) 

and at the interface, N] = -CD (dx;/dz) / (1 - xj) (3.52) 
the definition of h,, is —‘Ji* = Ay, (C1 i - Cha) (3.53) 
then hm C (%j- Xa) = -CD (dxj/dz) 

and dxj/dz = -h,, (xj- Xa) / D (3.54) 
finally Ny = C hy (4 - Xa) /(1 - x) (3.55) 
the Chilton-Colburn analogy is h / hy = Pair Cpair Le2/3 (3.56) 


and C is approximately constant (ideal gas behavior of air-vapor 
mixture) at the ambient value: C-=Caira=Pair/ Mair (3.57) 
also evap = Ni A M) and M) jj Mair = 0.624 (3.58) 


therefore evap = 0.624 (h A Le2/3 / cp air) (xj - Xa) /(1 - xi) (3.59) 


This detailed development reveals the sources of uncertainty in the 
evaporative term: the insolubility of air in the droplet and the Chilton- 
Colburn analogy and its 2/3 exponent. The property values of air, Cp air and 
Le, should be evaluated at the film temperature. The Chilton-Colburn 
analogy, derived by equating the functional forms of the equations relating 


Nusselt to Prandtl number and Sherwood to Schmidt number, makes the very 
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reasonable assumption that the small evaporative mass transfer rates do not 
influence the experimentally measured value of h. The liquid-vapor 


boundary condition can now be written 


-k 9T;/0z = h (Tj - Ta) + 0.624 (h A Le2/3 / cp air) (xi(T}) - Xa) / (1 - xi(T})) 
-F, 6 Tp4 [1 - 1/ (m(TR)d + D(TR) + 1] + He (3.60) 


The temperature dependence of xj; in (3.60) is evaluated by making the 
following assumptions: a) the interface is at thermodynamic equilibrium; b) 
the air does not alter the partial pressure of the water vapor at the interface; 
and c) the air and water form an ideal gas mixture. These approximations are 
very reasonable for the laboratory conditions and are frequently made in 


psychrometrics. Then 
xi(Tj) = Psat(Tj) / Pa (3.61) 


where Psat(Tj) is given for every degree Celsius in the Keenan steam tables 


[23]. 


Linearization 

Because (3.60) is nonlinear in T;, it cannot be directly used in the tridiagonal 
solution. Furthermore, (3.60) cannot be used to develop an analytical 
expression for the quasi-steady state temperature distribution. However, the 
code time step is chosen to allow only small changes in the liquid-vapor 
temperature between iterations. Then the temperature behavior of (3.60) can 


be represented by a line tangent at the current value of T;. This line is simply 


given by the first two terms of a Taylor series expansion: 
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Gcond(Tj) = qcond(Tjg) + (Oqcond(Tig)/ Tj) (Tj - Tio) (3.62) 
but Qcond(T}) = -k 0T;/dz = -k (AT; + B) = -k A Ls a k B 
and = qcond(Tj) = (Aqcond(Tjo) /ATj) Tj + Gcond(Tio) - (Oqcond (Tig) / Tj) Tio 


therefore A = -(Aqcond(Tjo)/0T;) / k < 0 (3.63) 
B = [(Aqcond(Tjo9)/2T;) Tip - qcond(Tig)] / k (3.64) 
where Acond(Tjg) = h (Tig - Ta) 


+ 0.624 (h A Le2/3 / cp air) (xi(Tig) - Xa) / (1 - xi(Tig)) 
-F, 6 Tr4[1 - 1 / (m(Tp)d + b(TR) + 1)] + He (3.65) 


The derivative may be found numerically, but the Clausius-Clapeyron 


equation provides an eloquent and reliable analytic expression: 


OQcond(Tj)/9T; = h + 0.624 (h A Le2/3 / cp air) AL (xi - xa) / (1 - xi)] /0T; 
but  ol[(xj- xa) / (1 - xj)J/0T; = {0[(xj - xa) / (1 - xi) /0x;} {0x;/0T;} 

= [(1- xa) / (1 - xi)2] Olpsat(T}) / pal /0T; 

= [(1 - xa) / pa(1 - xj)2] Opsat/OT; 

= [(1 - xa) / pa(] - xi)2] A / vig T; 
Oqcond(Tj)/OT; = h + 0.624 (h A2 Le-2/3 / pa cpair Vég Tj) (1 - xa) / (1 - x)? 


Therefore, the needed expression is 


Adcond(Tjg) /OT; (3.66) 
=h + 0.622 (h A? Le?/3 / ps Cp,air Vfg Tig) (1 - xa) / (1 - xi(Tig))? 


where, of course, Tjg must be in absolute units. 


As shown in Figure 24, the effect of the linearized boundary condition is to 
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FIGURE 24 
Liquid-Vapor Interfacial Flux vs. Liquid-Vapor Temperature 
(assuming 6 = 0.1 mm in the radiation term) 
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slide along the liquid-vapor flux curve and to overestimate slightly the 
updated T;. For T; greater than 95°C the curve becomes very steep and the 
linearization less accurate. However, the flux becomes rather stable at this 
point in the evaporative process, and only small temperature changes are 
expected. If the linearized boundary condition tries to assign an interfacial 
temperature greater than 100°C, then Tj is assigned a value of 99.9°C, and the 
liquid problem (tridiagonal transient or quasi-steady state) is solved using 


two boundary conditions of the first kind. 


Volume Flux 

Because the initial volume of the droplet Vg cannot be measured with great 
accuracy, there is little benefit to considering the small volume changes of the 
droplet associated with temperature changes. Therefore, the droplet density 
is assumed to be constant at the value of water at 20°C, 998 kg/ m3. Thus the 
droplet inventory can be made by volume rather than mass. In addition, the 
small effect of temperature on the axial grid size is also neglected. The 
following derivation gives the volume change calculated by EVAP during 
each iteration in terms of the previously derived evaporative flux (3.59). The 


summation is made over the NUMCOL annular columns of liquid. 


N1j = qevap,j / A Mi 

N1j Mi = p1 (AVj/At) / Aj 

AVj = devap,j At Aj / A p) 

but Aj=x { [j/NUMCOL]? - [(j - 1)/NUMCOL]? } Ro? 
= MR? (2j -1) / NUMCOL2 


AV = (nRg? At / A p} NUMCOL?) Sear (Qi) lasean (3.67) 


Again a special case must be made for the event of a false calculation of Ti 
greater than 100°C. The evaporative flux to be used in (3.67) is found from the 


known conductive, convective, and radiative fluxes: 


Gevap = qcond + rad - Gconv (3.68) 
where 

tridiagonal: qcond = -k (Tn - Tn-1) / Az (3.69) 
quasi-SS:_ — qcond = Hc + K(Tsg - 100°C - H.84/2k)/5 (3.70) 
and Gconv = h (100°C - Ta) (3.71) 
and Grad = Fe 6 Tr4 [1 - 1/ (m(Tp)d + b(TR) + 1)] - He 8 (3.48) 


3.8 Liquid Model as a Function of Radial Position 

The one-dimensional heat flux assumption limits the influence of one droplet 
on another to shape change by volume flux and to communication through 
the solid. The droplet shape (vector of annular droplet heights) has a strong 
effect on the time required until dT/ot can be neglected and the quasi-steady 
state solution substituted for the tridiagonal solution. This required time is 
quantified through the nondimensional time or Fourier number Fo = at/&2. 
The annular columns will become quasi-steady state starting with the 
thinnest, outermost column and proceeding inward. This is the rationale for 
the liquid solution structure as illustrated in the code flowchart (Figure 22). 
Provisions to change back to the transient solution are not needed because of 
the smoothly and slowly changing temperature profiles keep dT/ot 
negligible. The flowchart does have the provision for when the droplet is 
completely evaporated (or for geometry Model B, the droplet recedes past 
the column) that the liquid solution is bypassed and the solid-vapor 


boundary condition invoked. How to decide when to invoke the quasi- 
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steady state solution and the scaling problem of the outermost columns are 


discussed below. 


Inception of Quasi-Steady State Solution 


For the assumption of negligible residual volumetric heat generation, the 
linearity of the temperature profile indicates how close the column is to the 
quasi-steady state solution. The ratio of solid and vapor boundary fluxes 
indicates the linearity. Furthermore, the solid boundary flux starts high and 
decreases while the vapor boundary flux starts negative (radiation 
overpowers evaporative flux for early, low liquid-vapor interfacial 
temperatures) and approaches the solid flux. Thus the mechanism for 


changing the solution mode is 


FLUXR = (Tp - Tn-1) / (T2 - T}) 
IF (FLUXR .GT. PERC/100.) THEN (PERC set by user (0.97)) 


SSTi = (Tn + Tn-1) / 2 (transfer liquid-vapor temperature) 
NTRANC = NTRANC - 1. (decrement no. transient columns) 
END IF 


If He is not set to zero, then another approach must be used. Performing an 
investigative quasi-steady state solution on the outermost transient column 
during each time step and comparing solid boundary fluxes between the two 
models may be the most straightforward technique although directly 
calculating dT/dt values is a possibility. This point must be addressed if the 


code is to be used with nonzero values of He. 


Difficulty with Outermost Columns 


In order to capture the transient behavior of the thinnest, outermost columns, 
very small and computationally expensive time steps (order of 0.01 s) must be 
used; the Fourier number varies with the inverse square of 5. The BEM is 
extremely sensitive to time step size because it requires integration back over 
the previous time steps for each iteration. (The number of BEM computations 
varies with the inverse square of time interval, and the BEM accounts for 
more than half of the run time.) In fact, the outermost column may reach 
quasi-steady sate status before the heat wave reaches the droplet apex. If too 
large time steps are used in the Crank-Nicholson method, large erroneous and 
oscillating heat fluxes will be fed to the BEM destroying the credibility of the 
results. Thin droplets with large Bo are particularly a problem. Fortunately, 
because the transient behavior of these very thin columns is short, the BEM 


can tolerate flux estimates from these thin, transient columns. 


These thin columns are designated as quasi-steady state in the variable 
initialization stage of code execution. The Fourier number is the tool used to 
decide how thin a droplet must be to be assigned this special treatment. As 
the applicable nondimensional time, the Fourier number describes how far 
along in the transient process is the initially uniform temperature liquid 
column. A related problem suggests a threshold for Fo. Consider a column 
of unit thickness and unit diffusivity initially at temperature zero when one 
boundary is forced to temperature 1. Omitting the details, the transient 
temperature distribution can very accurately be represented by the first term 


of the eigenvalue solution: 
T(z, t) = 1-z-(2 /n) exp(-n? t) 
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The time until the midpoint reaches to within 98% of its steady-state value is 


calculated as 


1 - 0.5 - (2/n) exp(-n? t*) = 0.49 
t* = -In(0.01n/2) / n2 = 0.421 (3.72) 


The quantity t* is the Fo number for this system. Experimentation with 0.421 
as the threshold worked very well as there were no erroneous oscillations, and 
usually only the outermost one or two columns was affected. This solution to 
the thin-column problem works for any reasonable initial droplet geometry 


(Bp and Vo). 


A good initial liquid-vapor temperature is also required by the code to avoid 
invalid results. A poor initial value can result in the linearized boundary 
condition pushing the interfacial temperature past 100°C on the first time 
step, and recovery back down the steep flux curve (Figure 24) requires an 
unacceptable time. The following code is used to assign a sufficiently good 


initial liquid-vapor temperature: 


10 SST; = SST; + 0.1 
IF (LVFLUX(SST)}) - (Ts - SSTj)/8 .GE. 0.0) GOTO 20 
GOTO 10 
20 CONTINUE 


Again the case for nonzero He would add a complication here. 


3.9 Droplet Geometry - Models A and B 

As mentioned, the transient droplet shape must rely on an empirical model. 
Model A assumes a sphere segment shape and is defined by two inputs, the 
current volume V and the initial shape factor Bo, (or equivalently the wetted 


region radius Ro): 


6 = [current radius wetted region] / [radius of sphere with equal V] 


=2R/(6V/n)1/3 (defined by Bonacina [2]) (3.73) 


Because droplet surface tension is a function of temperature, and radiation is 
mostly absorbed at the droplet surface, the radiation can be expected to have 
a significant effect on Bp. However, EVAP takes Bo as a measured quantity 
and does not directly address this issue. Model B assumes a pancake shape 
and also allows input of the initial solid-liquid contact angle 69 and a | 


receding contact angle 0@;. 


Model A 

Photography of evaporating droplets suggests that the segment of sphere 
geometry may adequately describe the transient shape. Using analytic 
geometry the center of the sphere is adjusted axially until the sphere segment 
has the current volume V and the sphere boundary contains the solid-liquid- 
vapor interface. The expressions for 6(r) were developed by diMarzo and 


Trehan [5], checked by the author, and are repeated in the following: 


Y = (a/Ro) = [4/B3 + (1 + 16/B6)1/2J1/3 + [4/B3 - (1 + 16/B6)1/2]1/3 (3.74) 
5(r)/Ro = [((1/y + y)2/4 - (r/Ro)2]!/2 - (1/y- y)/2 (3.75) 
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The dependence of the contact angle 6 on the current value of B [5] is given in 
Figure 25. Model A updates the model shape when SUBROUTINE UPGEOA 
is called. 


Model B 

Evidence compiled by Chandra [24] shows that the droplet shape is somewhat 
flat, especially early in the evaporative process. Late in the evaporative 
process the wetted region begins to shrink as the receding contact angle 0, is 
reached, and the droplet surface tension and surface adhesion forces come 
into balance. EVAP uses a 6; value of 7°, a value suggested by Simon and 
Hsu who used experimental photographic techniques [25]. Model B has these 
two properties by making the following assumptions (as illustrated in 


Figure 26): 


(1) | When 86; is reached, the droplet shape may be considered to be a 
segment of a sphere. This assumption is reasonable because the 
spherical shape minimizes surface area before shrinking of the wetted 
region can begin. At the inception of this shrinkage the droplet has a 
unique volume. The aspect ratio remains constant as the droplet 
recedes. 

(2) The initial volume is the sum of a flat cylinder and a solid of rotation 
with outer surface tangent to the cylinder top face, with an acute 
contact angle 69 with the semi-infinite solid surface, and of circular arc 
profile. 

(3) | The gradual transition from initial shape to receding shape is uniquely 
specified by stipulating that the intersection of the horizontal tangent to 


the surface at the apex and the tangent to the liquid at the solid-liquid- 
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FIGURE 26 
Geometry Model B Transient Configuration of the Droplet 
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vapor contact point is always on a straight line. This convention 
reasonably approximates the experimental observations. 

(4) | The droplet apex will always be at or less than its initial value. Surface 
tension cannot force the droplet to pop up against gravity during 


evaporation. 


These assumptions have the effect of putting maximum and minimum limits 
on the initial contact angle 69. No shape meeting the four assumptions can 
contain the initial volume of liquid for values of 69 less than the 89 value 
given by the spherical segment configuration (Model A). Assumption (4) sets 
a maximum Q09 for the configuration with initial apex equal to the apex at the 
inception of shrinkage of the wetted area. Another upper bound for @09 is set 
at 90° by Chandra [24] for evaporating and boiling droplets (nucleate boiling). 
A geometrical calculation shows that for Bo less than 2.79, the maximum 69 as 
given by assumption (4) will always be greater than 90°. If one neglects the 
small recoil effect observed by Chandra [24], then 90° may be a good value to 
use. On the other hand, the @9 can be chosen as the value that yields the 
evaporation times and temperature profiles most closely matching the 
experimental data. Perhaps the Tso and impact velocity dependence of 89 can 


be discerned this way. 


Model B is invoked by calling SUBROUTINE UPGEOB, which was written 
by S. Tinker as part of an undergraduate project. It was slightly modified by 
the author to interface properly with the main routine, and the author 
thoroughly tested the subroutine for proper operation. The detailed 
workings of the 280-line subroutine are not explained here, but the code 


basically chooses among several potential types of droplet shape based on the 
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initial and current values of certain volumes. 


3.10 Coupling of the Liquid and Solid Domains 

The liquid and solid domains interact through the transfer of fluxes from the 
liquid to the solid and of solid-liquid interfacial temperatures to the liquid. 
The BEM solid solution uses fluxes at the center of the time intervals and 
gives surface temperatures at the end of the time intervals. The liquid 
solutions take the BEM temperatures and move up in time to the end of the 


time intervals producing fluxes valid at the ends of the time intervals. 


The simple predict-correct method, illustrated in the code flowchart 

(Figure 22), solves this mismatch in times while improving the reliability of 
the numerical results. Rolldown of the forcing function matrix FRCFNC 
refers to bumping values back in time as the new value is input. Rolldown is 
not made in the update of the forcing function matrix following execution of 
the predictor loop because the next BEM call will average the predicted 
fluxes and last corrected fluxes (from the same time step) to get a value 
representative of the center of the time interval. The resulting BEM 
temperatures are the corrected values. Because only the integration over the 
first recollection interval is new to this BEM call, significant computational 
time (on the order of 40% of the run time without this feature) is saved by 
using the integration values minus the first recollection step values (vector 
PREDU) in the corrected BEM call (SUBROUTINE BEM2). After time is 
incremented, the liquid corrector loop uses the corrected temperatures and 
the liquid temperatures from the last corrected loop execution to calculate 
the corrected fluxes. (The calculated predictor loop liquid temperatures are 


discarded.) These fluxes, valid at the beginning of the time interval, are given 


63 


to SUBROUTINE BEM1 for calculation of predicted surface temperatures, 


which are then inputted to the predictor loop. The cycle continues. 


3.11 Calculation of Property Values as Functions of Temperature 

Rather than pick representative values for most properties, the author chose 
to use the available computational power to use curve fit polynomials and 
function calls for property value evaluation. Then the effect of using 
representative values can be investigated and the use of property functions 
discontinued if warranted. However, it is the BEM and its thousands of calls 
to the bessel function sublibrary that is by far the most computational intense 
section of code. Following are brief descriptions of the sources and curve fits 


of the property data. Plots and tables are given in Appendix C. 


Lewis Number Le 

The Lewis number, a/D, is central to the evaporative boundary condition 
because it allows the use of the experimentally measured convective heat 
transfer coefficient h to calculate the rate of mass transfer. The air value for 


was constructed using data from [26] while D, given by [27], is 


D = 0.2232 (T/273)" cm2/s 273 <T < 370K (3.76) 


The resulting temperature dependence for Le (Figure 48) is small, and a 
representative value of 0.845 at 60°C is chosen. This value compares favorably 


with Le number data from another source. 


Constant Pressure Specific Heat of Air cp air 


The Chilton-Colburn analogy requires this quantity. Because it varies only 


0.7% over the range of interest, a curve fit is not required. However, Figure 49 


gives a five-point, quadratic least-squares fit anyway. The data is from [26]. 


Water Vapor Interfacial Mole Fraction xj 
The temperature dependence of this quantity, Psat/Pa, is crucial to the liquid- 


vapor boundary condition, so Table 2 was constructed for every degree C 
from the Keenan steam tables [23]. Figure 50 shows the exponential-like 


behavior of the interfacial mole fraction. 


Ratio of Latent Heat to Specific Volume Change A/vfg 

This combinaton is from the Clausius-Clapeyron equation and is also 
constructed from the Keenan steam tables [23] in Table 2. Figure 51 gives the 
temperature behavior. A/vfg and xj are input to EVAP through DATA 


statements. 


Latent Heat of Vaporization of Water A 


The six term curve fit formula given by Bejan [28] and suggested for 
computer applications is used. The formula and its constants are listed in 


FUNCTION HFG. 


Thermal Conductivity of Water _k) 


Values for k every five degree C from [16] were used in the fifth order 
polynomial fit of Figure 52. The thermal conductivity varies 13% over the 
range of interest. Of course, the linear quasi-steady state model demands a 
single value of k, and it uses the value at the current average of the boundary 
temperatures. The temperature difference between droplet boundaries will 


be seen to be small for the quasi-steady state solution. The error function 
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semi-infinite solution also uses such an average value. 


Thermal Diffusivity of Water a 
The thermal diffusivity, k/pcp, was constructed from the data in [16] using 


Table 3. It varies 15% over the range of interest, and the fifth order 
polynomial fit is given in Figure 53. Again an average value is used in the 


semi-infinite solution, particular to the user's choice of input temperatures. 


Density of Water p) 


The curve fit for the density of water is not used as was discussed in 
subsection 3.7; the value at 20°C, 998 kg/ m3, is used. However, the 
temperature dependence of the density of water is included in Table 3 [16] for 
the calculation of thermal diffusivity. The density varies 4% over the range of 
interest. The fifth order polynomial fit (Figure 54) is very accurate, more 
convenient than the available equations of state, and may find other 


applications. 


66 


4. SECTION III - RESULTS AND COMPARISON WITH LAB EXPERIMENTS 
After the code was fully debugged and in its final form, 20 runs were made in 
order to make comparisons with the single-droplet experimental results of 
Dawson [13] and perform a preliminary parametric study. The specific 


conditions for the experimental validation are 


Vo =40 ul; = Tso = 130°C; Tr. = 510°C; Bo = 1.61 


for which t=245 on average 


The Bo value is subject to significant experimental error because of the 
difficulty in discerning the droplet edge in Dawson's setup: the edge was 
implied by the temperature line scan data. Therefore, Runs 1 through 9 are 
devoted to the variation of Bp. Runs 10 and 11 investigate the effect of initial 
volume Vo while Runs 12 through 15 varied the initial temperatures. Use of 
geometry Model B was saved for Runs 16 through 19. Finally Run 20 
investigated a limiting case, that of a cylindrical droplet shape with constant 
wetted area. The code inputs (Tso, TL, Tr, geometry model, Vo, Bo, and 80 
(geometry Model B only)) and the principal output, the evaporation time 1, 
are given in Table 1. Also given are the time required for all annular columns 
to reach quasi-steady state status and the temporal/spatial averaged solid- 
liquid heat flux. The qdisk data refer to the constant heat flux model discussed 
in subsections 4.6 and 4.7. The figures in this chapter and Appendix D 


identify which inputs the graphed output is for by reference to run number. 


Figures 27 and 28 clearly show how the droplet steadily decreases in size until 


it is completely vaporized. The Model A droplet maintains its spherical 
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Results of Theoretical Calculations 


TABLE 1 
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segment shape at all times while the Model B droplet displays the more 


complex behavior for which it was designed. 


Figure 29 shows the cumulative energy transfers across the droplet boundary 
that must in the end add up to the latent plus sensible energy change in the 
original droplet mass. The conduction component is initially high because of 
the large temperature gradient at the solid interface. The direct radiation 
component is a almost straight line because it is not a function of any droplet 
temperature (it is a slight function of droplet thickness). The percentage of 
energy transfer by radiation is mostly just a function of wetted area and 
evaporation time since the flux is nearly constant. The percentage is 


approximately 15% for the run illustrated. 


Figures 30 through 40 display the temporal and spatial behavior of the solid- 
liquid flux, BEM forcing function, liquid-vapor temperature, and solid 
surface temperature for one case each of geometry Models A and B. The heat 
fluxes are highest at the droplet edge simply because the droplet is thinnest 
there. The effect is more pronounced for Model A as expected because of the 
larger very thin region. The forcing function graphs are just the flux graphs 
inverted and translated, but they do show how small the forcing function 
becomes outside the wetted area. There the only contribution is the 
convection and re-radiation due to the temperature depression from initial 
conditions. The oscillating behavior of the flux and forcing function for early 
times near the droplet edge originates from the discreteness of the numerical 


grid and would diminish with decreases in the radial and time steps. 


The liquid-vapor temperature graphs show that the upper surface is hotter at 
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Spatial Profile of Solid Surface Flux: Run 5 (Geometry Model A) 


—_—, 


SOLID SURFACE HEAT FLUX (kW/m’2) 


0 0.2 0.4 0.6 0.8 1 hr 1.4 
FIGURE 31 RADIAL POSITION r/R 


Spatial Profile of Solid Surface Flux: Run 16 (Geometry Model b) 
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FIGURE 33 
Spatial Profile of BEM Forcing Function: Run 16 (Geometry Model B) 
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LIQUID-VAPOR TEMPERATURE (C) 
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Spatial Profile of Liquid-Vapor Temperature: Run 5 (Geometry Model A) 
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Spatial Profile of Liquid-Vapor Temperature: Run 16 (Geometry Model B) 
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SOLID SURFACE TEMPERATURE (C) 
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Spatial Profile of Solid Surface Temperature: Run 5 (Geometry Model A) 
(t < T) 
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Spatial Profile of Solid Surface Temperature: Run 5 (Geometry Model A) 
(t > t; recovery) 
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SOLID SURFACE TEMPERATURE (C) 
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Spatial Profile of Solid Surface Temperature: Run 16 (Geometry Model B) 
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Spatial Profile of Solid Surface Temperature: Run 16 (Geometry Model B) 
(trecession < t < 7; droplet recession) 
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Spatial Profile of Solid Surface Temperature: Run 16 (Geometry Model B) 
(t > t; recovery) 
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the thin edge and that as the entire droplet reaches quasi-steady state status, Tj 
becomes fairly uniform. The solid surface temperature graphs show that for 
the combinations of inputs used the influence of the droplet is practically 
restricted to a disk of radius 4Rp. Additionally, after a short initial period, the 
solid surface temperature under the droplet becomes fairly uniform. This 
behavior was also produced in the theoretical results of Liao [12] and is due 
to the high fluxes at the droplet edge providing extra cooling to that annular 
region of the solid. Therefore, there is a straightforward relationship between 
solid surface flux and temperature as would be expected in a linear system. 
Figure 39 shows that the solid surface temperature profile has time to shift 
along with the droplet recession. In fact, Figures 37 and 40 show that the 


surface temperature begins to quickly recover after the droplet has vanished. 


The transfer of liquid solution control from the transient to quasi-steady state 
solution is illustrated by Figure 41. Every time the quantity FLUXR reaches 
0.97 the outermost transient column becomes quasi-steady state. The inner 
annular columns reach this status nearly at the same time because the droplet 


thickness is somewhat uniform there. 


4.1 Effect of Shape Factor Bo 

Figure 42 illustrates the effect of wetted surface area for the 10 yl droplet. The 
upwardly concave form for the curve indicates that the evaporation time 
decreases that come with increased surface area (goes with square of Bg) are 
not as large for large Bg. One possibility for this effect is that the nearly- 
constant radiation component has less time to make its contribution to 


droplets with lower evaporation times. 
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FIGURE 42 
Evaporation Time vs. Geometry Factor Bo: Runs 1-9 


The figures of Appendix D.1 illustrate more of the model behavior. The 
volume graph shows the steadiness of the dropwise evaporative process. The 
volume rate of change curves show how the evaporation mass transfer rates 
increase with liquid-vapor temperature. The nonsmooth features around 5s 
for high Bo possibly indicate numerical inaccuracies caused by the switch to 
the longer time interval. However, for small Bo the switch is remarkably 
smooth. The droplet bottom surface spatial average flux graph demonstrates 
the temporal behavior of the conduction contribution. After the high initial 
values, the average flux becomes roughly constant. Furthermore, the 
behavior is almost independent of Bop, meaning that the temperature 
difference across the droplet and droplet thickness effects cancel each other 
out. Figures 60 through 68 show precisely how the droplet upper and lower 


surface temperatures approach each other as the droplet thins. 


The recollection memory curves (Figure 69) show that solids with more 
compact droplets (small Bo) require less recollection time than solids with 
spread out droplets to know their present state. At the point that the droplet 
vanishes memory increases all the way back to deposition because the strong 
and recent forcing functions are gone and no longer overpower relatively 
long past events. The exponentially decaying behavior of the Green's 
functions makes it surprising that the memory should be as great as shown. 
The reason is that the Green's functions tend to level off for recollection times 
from 10 to 30 s and that the early forcing functions are the strongest; the effects 


compensate and the memory remains long. 


4.2 Effect of Initial Volume Vo 


The wetted area (and Ro) was kept constant in this investigation in an attempt 
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to learn independently the effect of the uncertainty in the experimental value 
of Vo. The effect was to cut t from 32 to 28 s for a 1 pl decrease in initial 
volume. The relationship between Vo and t is nearly linear for small changes 
at least in this case. Appendix D.2 gives the graphical results for this small 
parametric study. Figure 71 shows that the thinnest droplet (smallest Vo) has 
only a slightly higher mass transfer rate (as would be required for the linear 


relationship between Vo and 7). 


4.3 Effect of Initial Temperatures 

This investigation clearly shows the strong sensitivity of the solid-droplet 
system to initial solid surface temperature. Raising Tso 20°C from 120 to 140°C 
cut the evaporation time by almost 40%. The system is far less sensitive to 
changes in the initial liquid temperature and radiant panel temperature. This 
can be explained in that the sensible energy change in the droplet and the 
radiant contribution are relatively small parts of the required energy. 
Quantitative outputs are found in Appendix D.3. As expected the droplet 
initially at 60°C takes less time to get up to the maximum mass transfer rate 


than its 20°C counterpart. 


4.4 Effect of Geometry Model 

Implementation of geometry Model B produced some very interesting 
results, found in Appendix D.4. Figures 82 and 83 show that Model B results 
in longer evaporation times not because of the larger value of 69 but because 
of the slow down in volume flux after droplet recession. In fact, Figure 84 
shows that time to reaching the receding contact angle is nearly independent 
of initial contact angle. The similar plot with time scaled to t (Figure 85) 


shows that the effect of increasing Bo is to reach recession comparatively 
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faster. The volume plot (Figure 82) suggests that the special cylinder- 
geometry-model droplet would have the longest evaporation time of the Bo = 
2.0 droplets if it had a receding property. This coupled with the observed 
order of the curves leads to the conclusion that evaporation time increases 
slightly with 80, all other quantities equal. Severe local thinning at the edge, 
where area is greatest, is more efficient at pushing flux through the wetted 


area than moderate thinning throughout the droplet. 


4.5 Comparison with Experimental Results 
Comparisons between the theoretical and experimental results can be made in 
two ways: through the evaporation time and through the radial and temporal 


solid surface temperature dependence. 


Evaporation Time Comparison 
If one accepts all the values given by Dawson [13], then geometry Model A 


predicts an evaporation time 33% too long and Model B 58% too long. 
However, there are several reasonable explanations that the source(s) of the 
discrepancy may be related to the experiment (possible sources of 


discrepancy related to the code are discussed in Chapter 5): 


(1) Kidder [8] performed experiments very similar to Dawson's and 
measured an average t of 31.2 s for an average Bg of 2.0. This is very 
close to the 32 s of Model A. Kidder's radiant panels were of a different 
geometry, but this effect can be argued not to be significant for several 
reasons: a) Chapter 2 showed that the reflectivity dependence on 
incident angle cannot be expected to change Fg much; b) overall energy 


balances and equal Tso values mean that the radiative fluxes should be 
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close; and c) the spectral dependence of Kidder's incident radiation (TR 
value) would not be too different from Dawson's. Then, the simplified 
radiation model predicts very similar droplet behavior. 

(2) Bo was difficult to experimentally measure in Dawson's setup. Also 
radiation-induced surface tension relaxation at the droplet edge may 
mean that Ro was really larger than it seemed. 

(3) Slight errors in the measurement of Vo or Tso could have significant 
effect on evaporation time. 

(4) | Dawson measured evaporation time by marking the time at which the 
video tape showed the temperature profile begin to snap back toward 
the initial conditions. Results using geometry Model B showed that the 
snapping back of the temperature profile would begin at the point of 
recession not complete evaporation. Figure 82 shows that this effect 


may be on the order of 4 or 5s. 


Solid Surface Temperature Profile Comparison 


This comparison is greatly complicated by the difference in evaporation 
times and uncertainty with Bo. The Bo = 2.0, Model A theoretical results were 
used for comparison because of the agreement with Kidder's value of Bo. 
Discrepancies in Bo values cause a scaling problem with the experimental 
data, which is given as a function of r/Ro. The purpose of Figures 43 and 44, 
which show the temperature profiles overlayed, is only to demonstrate the 


agreement of the general behavior of the theoretical and experimental Ts. 


4.6 The Constant Heat Flux Approximation 
The spatial and temporal behavior of the solid surface fluxes suggest that the 


classical constant heat flux solution might be a useful model of the 
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Run 5 Temperature Profile Compared with Experimental Data [13] (t < 7) 
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FIGURE 44 RADIAL POSITION r/R 
Run 5 Temperature Profile Compared with Experimental Data [13] 
(t > t; recovery) 


evaporative process. Carslaw and Jaeger [29] and later Beck, et al., [20] 
showed that for a semi-infinite solid surface suddenly subjected to a constant 
and uniform heat flux q out of solid over a circular region of radius R with no 


other heat transfer on the remainder of the surface: 


Tso - T(r, 0, t) = (qR/ks) J.” Jo(Br/R) Jx(B) (1/B) erf [B(ast)/2/R] dB (4.1) 


This expression can be developed using the Green's function solution 
equation (GFSE) and judicious analytical integration. The integration on B 
originates from the particular Green's function expression most amenable to 
the problem. This solution has no provision for the case when the droplet is 


evaporated and the surface temperature begins to recover. 


The author overcame this limitation by changing the limits on the recollection 
time integral of the GFSE from (0, t) to (t - t, t) where t - t is the time since the 
flux disk was suddenly removed. After lengthy manipulations the solution 


valid for times before and after disk removal is obtained: 


Tso - T(r, 0, t) = 
(qR/ks) Pe Jo(Br/R) Ji(B) (1/B) {erf [B(ast)1/2/R] - erf [B(astt)!/2/R] } dB 


where tt=0 fort <tT 


(tz (ee fort >T (4.2) 


Note that as t approaches infinity both error functions approach 1 and cancel 


each other so that the initial temperature Tso is restored. 
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4.7 Comparison with the Constant Heat Flux Model 

Because there is an initial flux in the droplet-solid system, the value q in (4.2) 
must be defined to include the effect of subtracting out the initial conditions. 
Recall that the full GFSE includes a term for initial conditions that can be 


ignored only by mapping to another domain. The required definition for q is 


q =disk- qo where qdisk > 0 and qo < 0 for the radiation case (4.3) 


The initial flux qo is obtained by an overall energy balance. The conductive 
flux qdisk is an output of EVAP that may be curve fit for conditions of interest 
to use in (4.3) and (4.2). Alternatively, the theoretical or experimental 
evaporation times can be curve fit, and an energy balance around the droplet 


used to find qdisk: 


disk ™ R2 t + Grad T R2 t - Gcony 7 R2 t = pj Vo Ah 
disk = P1 Vo Ah / x R21 + Gconv ~- rad (4.4) 
where Anz Cp AT+A 


and AT is a representative value of the droplet temperature change before 


evaporation. 


The computer code needed to implement (4.2) is given in Appendix E. The 
subroutine QAGI is called to a system library in order to perform the difficult 
semi-infinite numerical integration. The argument is a bessel-type damped 
oscillation which is difficult to integrate conventionally. In addition, the 


subroutine requires the code to be compiled in double precision mode. 


Typical results are given in Figures 45 and 46. The curves compare very 
favorably with the EVAP theoretical and experimental data. Note that the 
temperature profile plateaus under the droplet are curved for the constant 
flux model. This is expected because the high edge fluxes are brought down 
to the average. The inability of the model to duplicate the behavior of droplet 
recession, especially for high Bo (early recession), is a weakness. However, an 
additional curve fit to theoretical or experimental data by mapping t to some 
function of t in (4.2), i.e. replacing t with f(t), could improve the model 
validity. Another approach is to superimpose two concentric disks of 
variable relative strengths or areas. Use of a validated constant flux model is 


the current approach for a multi-droplet code in this line of research. 


4.8 Use of the Constant Heat Flux Model to Test the BEM 

As was mentioned in the subsection on the BEM, the constant heat flux model 
was used to validate the BEM section of EVAP. The procedure was simply to 
redefine the forcing function to be 1 forr <1and0forr>1. All other 
quantities were also nondimensionalized. The nondimensional form of (4.2) 
was also computed. The nondimensional temperature profiles are compared 
at three times in Figure 47. The maximum deviation is about 1.5% of the full- 
scale temperature and occurs around anr value of 0.5. This test says nothing 
directly about how adequate the BEM grid size is for r > Ro because the 
forcing function is zero there. However, the forcing function is normally 
small for r > Rg, and the grid size is Ro/12 for r < 4Ro. Therefore, the author 


concludes that the BEM section of EVAP is functionally well. 
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FIGURE 45 
Constant Heat Flux Model (t < 7) 
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Constant Heat Flux Model (t > t; recovery) 
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FIGURE 47 
BEM Section Validation Using Constant Heat Flux Model 
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5. CONCLUSIONS AND RECOMMENDATIONS 

The code output clearly demonstrates that the program is functionally 
smoothly with no numerical glitches. However, some further 
experimentation may help to shed light on the experimental deviations. 
Other combinations of Tso, Tr, Bo, and grid size will surely help. It is 
improbable that the radiation model is the source of the deviation because of 
its relatively minor role for the Tso = 130°C case; increasing the radiative flux 
50% in Run 13 did not change the evaporation time. The geometry models are 
always questionable because they are empirical simplifications. Attempting 
to validate the program using Kidder's data [8] may be the most 
straightforward approach. The code was written to ensure that future 


researchers could carry on the work without a steep learning curve. 


Although this effect did not appear in the results given in Chapter 4, the 
EVAP code has been observed to grossly overestimate the wetted region heat 
flux for the last time step with evaporation. This is possible if the droplet 
vanishes in the early portion of this time step. Because the droplet is very thin 
in this case, the calculated flux is very high and not representative of the 
entire time step. The remedy is to replace the false high flux with the uniform 
flux that would be needed to evaporate the remaining liquid over the entire 
time step. The following should be inserted immediately above the updating 


of the droplet volume: 
IF (V .LT. DV) THEN 
DO 185 J = 1,.NUMCOL 


IF (DEL(J) .GT. 0.) 
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+ FLUX(J) = V * RHOH20 * HFG(SSTI(J)) / (PI * RADIUS**2 * DT) 
185 CONTINUE 
END IF 


The expression for the constant heat flux model helps to explain one 
difference between the cases of heating by conduction from below and by 
radiation from above. For equal but opposite initial fluxes qo, the q quantity 
of (4.3) has a value 2qo larger for the radiation case. Thus the temperature 
depression is greater for the radiation case, and the temperature at the solid- 
liquid interface can be expected to be lower. Another point of view is that the 
heat flux has the more difficult path to go from the solid surface down, 
radially in, and then up into the droplet. In the conductive case the path is 
straight up into the droplet. The effect of direct radiation to the droplet and 
radiation-induced surface tension relaxation of the droplet (increased f) 
counteract this enhanced temperature depression effect. The result of these 


compensating effects is rather similar evaporation times. 


Future uses for the code with minor modifications include the following: 

(1) Calculate the H(z) and F(z) functions for Tr values beyond those easily 
possible in the laboratory using the radiation model code. Pick 
nonzero Hc values and investigate the effect of high temperature 
incident radiation. The magnitude of the radiative flux can be kept 
reasonable by using small view factors. 

(2) Radically alter the liquid material properties in order to test the 
cooling behavior of other liquids or hypothetical materials. 

(3) Investigate in detail the effect of the solid properties (thermal 


conductivity and diffusivity). 
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Appendix A QuickBASIC Code for Radiation Model 
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REM Program FRAD.BAS calculates the Radiative Flux Frad 
REM as a function of z (distance to surface) given the file FTHETA.DAT 


REM for ftheta using the trapezoidal rule of double integration 


CLS 

DEFDBL A-H, O-Z 

CONST pi = 3.14159265358979# ; 

DIM NumLam(1 TO 3) ‘Highest Number for Jth Segment of Integral 


DIM dlam(1 TO 3) ‘wavelength step size in um for Jth Segment 


numloop = 0 


INPUT “Enter the Coil Temperature Tcoil (C): J Gitexe lil 

INPUT “Enter the local distance from the surface z (mm): a) We: 
Teoil = "Tcoil + 273.15# 

Gnhe= 2334 

Cl = 374200000# 

C25= 143:90# 


REM Read in ftheta values 
DIM fdata(0 TO 900), thetadata(0 TO 900) 
OPEN "A:\FTHETA.DAT" FOR INPUT AS #33 
FOR K = 0 TO 900 

INPUT #33, thetadata(K), fdata(K) 
NEXT K 
CLOSE #33 


dtheta = pi * .1# / 180# 
NumLam(1) = 16: NumLam(2) = 20: NumLam(3) = 10 
dlam(1) = .05#: dlam(2) = .2#: dlam(3) = .5# 


SUM3 = 0# 
FOR J = 1 TO 3 ‘Three Integration Segments are Necessary 
SUM2 = 0# 
FOR Ilam = 0 TO NumLam(J) 
SUM = 0# 
READ wl, alam 
tau-earam * (Zz / 1LO#') ‘z in mm, alam in cm*-1 
Pea Ce mew ot. = (EXP(C2) /e (wl, * Tcorl))e— 1)) 
PRINTS oeSe is, Wl =" wl, “taue= “> tau 


SUM = SUM + O# ‘First Value is 0 
FOR thetad = .1# TO 89.95# STEP .1# 
theta = pi * thetad / 180# 


REM Retrieve ftheta 

K = INT(thetad * 10 + .5) 
ftheta = fdata(K) 
thetainput = thetadata(K) 


REM Calculate rho 

US=N (Chane 2#e—eSIN (theta) a) 2#)) S55 ¢ 

Rperp = ((cOS(theta) - u) / (COS(theta) + u)) * 2# 

Rpame=(Cneoez#. 2s COS (theta) — Ww) / (Cn o 24) *eCos( thetaye+ nu) > 2¢# 
rho = (Rperp + Rpar) / 2# 


REM Calculate mu (cmwu here) 
A = .75# * SIN(theta) 
cmun=—COS(ATIN(A/) (a 5-0 Ay) 24) = 2 5#)) 


PRINT thetad, thetainput, ftheta, K, rho, cmu 


funct = ftheta * COS(theta) * SIN(theta) * (1# - rho) * EXP(-tau / cmu) 
SUM = SUM + 2# * funct 
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SUM = SUM + O# Pat: meus pas °, : Y 
: : ‘E not included in func 
Sn SUM ap om hom Lettatts ales ene 


PRINT "SUM = "; SUM, "SUM2 = "; SUM2, "SUM3 = "; SUM3 


, 
PRINT TrrTrrcrereee eee ee eee eee eee 


FOR Idelay = 1 TO numloop: NEXT Idelay 
IF Ilam > O AND Ilam < NumLam(J) THEN 
SUM2 = SUM2 + SUM * dtheta 
ELSE 
SUM2 = SUM2 + SUM * dtheta / 2# 
END IF 
NEXT Ilam 
PRINT “ttt e$HHHHHHHH HHH HEHEHE HEHEHE HEHEHE ++ ++ Ht++” 
PRINT "SUM2 = “; SUM2 
PRINT "HttHEHEEEHEEEHE HEHEHE HEHEHE HEHEHE HEHEHE HH HHH H++" 
FOR Idelay = 1 TO 2 * numloop: NEXT Idelay 
SUM3 = SUM3 + SUM2 * dlam(J) / 2# 
NEXT J 


PRINT *SUM3 = *; SUM3 
PRINT 


Frad = 2# * SUM3 / 100# 

PRINT “Frad = "; Frad; “ kW/m*2" 
PREN TS Ze = eee eee tN ue 

s}sGqipi 4rvefepiil <3 Wo yejeslale GY (eh 


REM Absorption coefficient of water alam (cm*-1) as function of wavelength 
DATA) .20, %. 0697 
DATA .25, .0168 
DATA p20 npr 00G7 
DATA .35, .0023 
DATA .40, .00058 
DATA .45, .00029 
DATA .50, .00025 
DATA .55, .000045 
DATA .60, .0023 
DATA .65, .0032 
DATA .70, .0060 
DATA. a7 Shey. 0261 
DATA .80, .0196 
DATA .85, .0433 
DATA .90, .0679 
DATA .95, .388 


DATA 1.0, .363 
DATA 350; .363 
DATA 1.2, 1.04 
DATA 1.4, WAL 
DATA 1.6, 6-72 
DATA 1.8, 8.03 
DATA 2.0, 69.1 
DATA 2.2, 16.5 
DATA 2.4, 2K0) al 
DATA 2.6, 153 
DATA 2.8, 5160 
DATA 3.0, 11400 
DATA 3.2, 3630 
DATA 3.4, nb 
DATA 3.6, 180 
DATA 3.8, Tre 
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lculates the Volumetric Heat Generation Rate Hrad 
(distance to surface) given the file FTHETA.DAT 
he trapezoidal rule of double integration 


REM Program HRAD.BAS ca 
REM as a function of z 
REM for ftheta using t 
see H, O-Z 

LwA=H; - 
Sone: pl = 3.14159265358979# 
DIM NumLam(1 TO 3) 
DIM dlam(1 TO 3) 


‘Highest Number for Jth Segment of Integral 
‘wavelength step size in um for Jth Segment 


numloop = 0 


INPUT “Enter the Coil Temperature Tcoil (C): Os Uterow ish 

INPUT "Enter the local distance from the surface z (mm): 7 ae a 
Tcoil = Tcoil + 273.15# 

Cnig=eleoed 

C1 = 374200000# 

C2=514390F 


REM Read in ftheta values 
DIM fdata(0 TO 900), thetadata(0 TO 900) 
OPEN "A:\FTHETA.DAT" FOR INPUT AS #33 
FOR K = 0 TO 900 

INPUT #33, thetadata(K), fdata(K) 
NEXT K 
CLOSE #33 


dtheta = pi * .1# / 180# 
NumLam(1) = 16: NumLam(2) = 20: NumbLam(3) = 10 
dlam(a)ia= )0S#:ed Lam (2) = se2h cdl am (3) =e or 


SUM3 = O# 
EORw: =p ‘TOms ‘Three Integration Segments are Necessary 
SUM2 = 0# 
FOR Ilam = 0 TO NumLam(J) 
SUM = 0O# 
READ wl, alam 
talie= alames (29/7) LO) ‘z in mm, alam in cm*-1 
BE =sCl (awh Ge5#  s. (EXP C270/ (wile TCO 1) )) emu), ) 
PRING. “Ee=e20 <b, OWle> “wiles otauec! cr etau 


SUM = SUM + 0# ‘First Value is 0 
FOR thetad = .1# TO 89.95# STEP .1# 
theta = pi * thetad / 180# 


REM Retrieve ftheta 

K = INT(thetad * 10 + .5) 
ftheta = fdata(K) 
thetainput = thetadata(K) 


REM Calculate rho 

u = (Cn * 2# - SIN(theta) * 2#) * .5# 

Rperp = ((COS(theta) - u): / (COS(theta) + u)) * 2# 

Rpar’= ((Cn ee2# * COS( theta) ~-91) 7 7 (Che 28 GCOS (Chetay erat) ee eae 
rho = (Rperp + Rpar) / 2# 


REM Calculate mu (cmu here) 
A = .75# * SIN(theta) 
cm, =" COS(ATN(A 7 (1 =9A 1424) 5 #)) 


PRINT thetad, thetainput, ftheta, K, rho, cmu 
funct = ftheta * COS(theta) * SIN(theta) * (1# - rho) * EXP(-tau / cmu) 
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SUM = SUM + 2# * funct 
NEXT thetad 


SUM = SUM + O# ‘Last Value is 0 

SUM = SUM * E * alam "E not inctuded in funce 
PRINT NEA ‘ 

PRINT I SUMe= ae SUM, s USUM2a=9 >) SUM2, 9SUM3 = si SUMS 


PRINT eee eee eee eee eee 


FOR Idelay = 1 TO numloop: NEXT Idelay 
IF Ilam > 0 AND Ilam < NumLam(J) THEN 
SUM2 = SUM2 + SUM * dtheta 
ELSE 
SUM2 = SUM2 + SUM * dtheta / 2# 
END IF 
NEXT Ilam 
PRINT "HEHEHE HEHEHE HEHEHE HEHEHE HHH HH HH HH +++ te+e +" 
PRINT "“SUM2 = "; SUM2 
PRINT “He tete He tHeee HHH ttt te Hee eee eee eetetett+++" 
FOR Idelay = 1 TO 2 * numloop: NEXT Idelay 
SUM3 = SUM3 + SUM2 * dlam(J) / 2# 
NEXT J 


PRINT “SUM3 = "; SUM3 
PRINT 


Hrad = 2# * SUM3 

PRINT “Hrad = "; Hrad; " kW/m*3° 
PRUN Dee cs iemee cee) Is 

PRIN CO m=—mrorreT COL ss © (c" 


REM Absorption coefficient of water aiam (cm*-1) as function of wavelength 
DATAG2 2070070691 
DATA .25, .O0168 

DATAT 3.0 ee O0I6i7 

DATAS 35,8 Woes 

DATA .40, .00058 

DATA .45, .00029 

DATA .50, .00025 

DATA .55, .000045 

DATA) 200),...0023 

DATA .65, .0032 

DATA .70, .0060 

DATA fa, 0261 

DATA .80, .0196 

DATA, 485;) .0433 

DATA .90, .0679 

DATA .95, .388 

DATA 1.0, .363 


DAT ANS 08 .363 
VATA 2.2, 1.04 
DATA 2.4, 12.4 
DATA 2.6, 6.72 
DATA 1.8, 8.03 
DATA 2.90, 69.1 
DATA 232), 16.5 
DATA 2.4, = Ore 
DATA 2.6, iB. 
DATA 2.8, 5160 
DATA 3.0, 11400 
DATA 3.2, 3630 
DATA 3.4, 721 
DATA .3.6;, 180 
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DATA 
DATA 
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DATA 
DATA 
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Appendix B FORTRAN Code for Single-Droplet Model 
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3 SLENN WHITE 

¢c (301) 405-5334 

fad ADVISOR: MARINO DIMARZO (301) 405-5257 

Cc MASTER'S THESIS -- UNIVERSITY OF MARYLAND AT COLLEGE PARK 

c EVAPORATIVE COOLING WITH RADIANT HEAT INPUT 

(e PROGRAM TO CALCULATE TRANSIENT SURFACE TEMPERATURE DISTRIBUTION 
c 

c 

C2 


WRITTEN FOR NIST CENTER FOR FIRE RESEARCH 


eT: DEFINE COMMON BLOCK TO ALLOW SUBROUTINE UPGEOB TO KEEP VALUES OF CERTAIN VARIABLES STATIC 
REAL BVAR1, BVAR2, BVAR3, BVAR4, BVARS, BVAR6, BVAR7, BVARS, BVAR9, BVAR1O 
REAL BVAR11, BVAR12, BVAR13, BVAR14, BVAR1S 
COMMON /GEOBBL/ BVAR1, BVAR2, S3VAR3, BVAR4, BVARS, BVAR6, BVAR7, BVAR8, 
BVAR9, BVAR1O, BVAR11, BVAR12, BVAR13, BVAR14, BVAR1S 
PARAMETERS WHICH MAY BE SET HERE BUT ARE NOW INPUT AS VARIABLES THROUGH THE MONITOR.......... 


Some 0 

Ch avetey s THE TIME TO EXIT THE TIME STEP LOOP TEND IN S 

c PARAMETER (TEND = 60.) 

Gnecxshe THE INITIAL VOLUME VO IN M**3 

c PARAMETER (VO = 10E-09) 

Gorter THE INITIAL SHAPE PARAMETER BETAO 

c PARAMETER (BETAO = 2.00) 

Grteteier THE GEOMETRY MODEL DESIRED, ‘A’ OR ‘B’ 

c PARAMETER (GEOMOD = 'B’) 

Cre THE INITIAL CONTACT ANGLE THETAO IN DEGREES (>THETAA) (GEOMOD = ‘B’ ONLY) 

i PARAMETER (THETOD = 60.) 

(Shere ete THE INITIAL SURFACE TEMPERATURE (TSO) 

Cc PARAMETER (TSO = 130.) 

Starck THE RADIATIVE HEATER COIL TEMPERATURE IN C (475-650) 

Cc PARAMETER (TCOIL = 510.) 

Givers NEGLECT RATIO TO DECIDE WHEN TO TERMINATE MARCH BACK INTO TIME 

c PARAMETER (NEGRAT = 0.01) 

Creer ROOT OF OUTPUT FILE NAMES (6 CHARS) 

Soren FILE NAME OF OUTPUT OF MAIN DROPLET PARAMETERS IS (FILNAM//' .dat’) 

(Shige Ac ILE NAME OF OUTPUT OF SOLID SURFACE TEMPERATURES AND FLUXES IS (FILNAM//' .out’) 
(SS, FILE NAME OF OUTPUT OF SPATIAL DROPLET SURFACE TEMPERATURES AND FLUXES IS (FILNAM//’ .ins’) 
Son FILE NAME OF OUTPUT OF SPATIAL DROPLET HEIGHT AS FUNCTION OF TIME IS (FILNAM//' .geo’) 
Car FI NAME OF OUTPUT OF CUMULATIVE ENERGY TRANSFERS IS (FILNAM//’' .bal’) 

cS PARAMETER (FILNAM = ‘2run0l’) 

Ss NUMBER OF TIMES TO OUTPUT SPATIAL TEMPERATURE AND FLUX DATA 

= PARAMETER (NDATIM = S) 

C.....SPECIFY TYPES OF THE INPUT VARIABLES DEFINED ABOVE 


REAL TEND, VC, BETAO, TSO, TCOIL. NEGRAT, THETOD 
CHARACTER*6 FILNAM 

CHARACTER®10 FILNA1, FILNA2, FILNA}. FILNA4. FILNAS 
CHARACTER*®1 GEOMOD 

INTEGER NDATIM 


PARAMETERS WHICH CAN BE ADJUSTED HERE WITHOUT ANY CTHER CHANGES ---------------------- 522 rr nnn -- 
THE RECEDING ANGLE THETAR IN DEGREES GEOMOD - B’ ONLY) 

PARAMETER (THETRD = 7 } 

THE AMBIENT TEMPERATURE 

PARAMETER ‘TAMB = 25 ) 

THE INITIAL DROPLET TEMPERATURE 

PARAMETER (TL = 25.) 

MINIMUM FLUX RATIO PERC TO DEFINE SS :NA [Ci ES 99) 

PARAMETER (PERC = 97 


te} a} an 


ia) 


PARAMETERS NOT TO BE CHANGED WIZ Sw STMER SHANGES << <2 2 26 sence eee cre mesene men sasesciecccscesccic= 
WARNING: DO NCT ADJUST TIMES EXCEPT TEND: WITHOCT CHANGING SUBROUTINES RECONF, WEIGHT, AND WGHT 
THE STARTING TIME STEP OTSHRT IN S 

PARAMETER (OTSHRT = 0.1) 

THE TIME STEP DOTLONG IN S 

PARAMETER (DTLONG = 1.0) 

THE TIME TO USE THE STARTING TI STEP TSHORT <N S 


ana 


fa) 


n 
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PARAMETER (TSHORT = 4.’ 
S NUMBER OF ANNULAR COLUMNS TO MCDEL THE DROPLET 


PARAMETER (NUMCOL = 12) 
TOTAL NUMBER OF NODES N PER COL (DZ=DELTA/ (N-2)) 


PARAMETER (N = 12) 
THE NUMBER OF TIMES USED FOR THE BEM SUBROUTINE 


C ARAMETER (NUMBET = 100) 
re HE NUMBER OF BEM NODES USED 
PARAMETER (NUMNOD = 78) 
ee PARAMETER (PI = 3.14159265358979) 
c E TYPE Vp qed aes bh doh NR EPO a ERI ied BAB dt AN PS ACME COTS Eh ROCK OCB CN 
REAL GAMMA, THETAO, THETAR. V, RADG, RADIUS, VR, RHOH20 
REAL ELSIS, ALFSIS, EPSF, SIGMA. F, EPSILN, FRADO, FRADS 
REAL HR. H, HO, FLUXO, TC. OT, KS, FLUXOS, TIME, TI 
REAL NFLUX, DVCOL, DV, RAV. QC, SSFI, SSDFDT, SSA, SSB, SSTEMP 
REAL EXPR, SSTOLD, FLUXR, SSTIME, EVTIME, THTIME, MEMORY, BETA, THETA, FO 
REAL CONDUC, RADIAT, CONVEC, LATENT, MAXSEN, SENSIB 
REAL TICONV, TICAVG, TAVG, FAVG 
INTEGER I, J, K, L, NTRANC 
CHARACTER*1 TSFLAG, TRNBGN, EVFLAG, THFLAG 
te DESCRIPTION OF VARIABLES (ALL UNITS ARE SI EXCEPT TEMPERATURE IN C): 
une GAMMA: CRANK-NICHOLSON CONSTANT; THETAO: INITIAL DROPLET CONTACT ANGLE; 
Ramet THETAR: RECEDING CONTACT ANGLE (GEOM MODEL 8 ONLY); V: CURRENT DROPLET VOLUME; 
Che RADO: INITIAL WETTED RADIUS; RADIUS: CURRENT WETTED RADIUS; 
ae VR: VOLUME AT WHICH RECEDING BEGINS (GEOM MODEL B ONLY); 
Paiey RHOH20: DENSITY OF WATER; KLSIS: SEMI-INFINITE SOLUTION (SIS) VALUE OF CONDUCTIVITY FOR DROPLET; 
aes ALFSIS: SIS THERMAL DIFFUSIVITY VALUE OF DROPLET; 
Rees EPSF: ABSORPTIVITY OF DROPLET GIVEN DISTRIBUTION OF INCIDENT RADIATION BY POLAR ANGLE; 
ae 8 SIGMA: RADIATION CONSTANT; F: THEORETICAL TILE VIEW FACTOR; 
“8 4G EPSILN: SOLID EMMISIVITY; FRADO: RADIATIVE FLUX PENETRATING INTO SURFACE OF DROPLET; 
ae FRADS: RADIATIVE FLUX ABSORBED BY SOLID; HR: RE-RADIATION HEAT COEFF OF SOLID; 
eet H: CONVECTIVE HEAT COEFF; HO: OVERALL HEAT COEFF; FLUXO: INITIAL HEAT FLUX; 
eer TC: SIS CONTACT TEMPERATURE; DT: TIME STEP; KS: SOLID CONDUCTIVITY; 
Ck: FLUXOS: SIS HEAT FLUX FOR TIME STEP 0.5; TIME: ELAPSED TIME SINCE DEPOSITION; 
ee TI: LIQUID-VAPOR INTERFACIAL TEMPERATURE; NFLUX: MOLAR FLUX LEAVING DROPLET; 
Ane DVCOL: CHANGE IN VOLUME OF ANNULAR DROPLET COLUMN; DV: CHANGE IN VOLUME OF DROPLET; 
ra KAVG: AVERAGE STEADY-STATE (SS) CONDUCTIVITY; QC: CONDUCTIVE HEAT FLUX; 
Canoes SSFI, SSDFDT, SSA, SSB, SSTEMP, EXPR, SSTOLD: SS PARAMETERS FOR TAYLOR LINEARIZATION OF LIQUID-VAPOR BC; 
eee FLUXR: RATIO OF BOUNDARY TEMPERATURE GRADIENTS OF OUTERMOST TRANSIENT ANNULAR DROPLET COLUMN; 
he SSTIME: TIME UNTIL ALL COLUMNS REACH SS; EVTIME: TIME UNTIL DROPLET EVAPORATES; 
eee THTIME: TIME UNTIL RECEDING ANGLE IS REACHED (GEOM MODEL B ONLY) 
ae MEMORY: RECOLLECTION TIME AT WHICH EFFECT OF FORCING FUNCTION BECOMES NEGLIGIBLE (BASED ON NEGRAT) ; 
on BETA: CURRENT SPLAT FACTOR (BASED ON RADIUS AND SPHERE OF CURRENT VOLUME); 
a ee THETA: CURRENT DROPLET CONTACT ANGLE IN RADIANS (GEOM MODELS A AND B); 
can. FO: FOURIER NUMBER ALPHA*TIME/DEL(J)“2 USED TO DETERMINE IF COLUMN SHOULD BE SS INITIALLY; 
oe CONDUC, RADIAT, CONVEC, LATENT. MAXSEN. SENSIB: CUMULATIVE ENERGIES; 
aA Se TICONV: LIQUID-VAPOR INTERFACIAL TEMPERATURE; TICAVG: AREA AVERAGED TICONV; 
pbs TAVG: AREA AVERAGED TEMPERATURE UNDER THE DROPLET; 
Gate FAVG: AREA AVERAGED HEAT FLUX UNDER THE DROPLET; ; 
ei I. J, K: ARRAY INDICES; NTRANC: CURRENT NUMBER OF TRANSIENT COLUMNS; 
C.....TSFLAG: FLAG -- ‘1° = (TIME > TSHORT), 
Bi cate TRNBGN: FLAG -- ‘1’ » THE TRANSIENT TRIDIAGONAL SOLUTION METHOD HAS BEGUN ALREADY; 
cae EVFLAG: FLAG -- ‘1’ = (TIME > EVTIME); 
eos THFLAG: FLAG -- ‘1’ » (TIME > THTIME) 
S$ a eens [RRO RRRRRRRRERRERRRRRRRER RE RR RRR RE RRR RR RRRER RR RRRRR RRR RRRR RRR RRR RRR R RR RRRARR RASS SS 
Coe FUNCTION TYPE SPECIFICATIONS 
REAL ALPHAL, CTHETA, CPA, DFDT. FRAD, HCONV, HFG, HV, KL 
REAL LVFLUX, NH20, SISTC, SOLVAP. XI 
c Pr DIMENSION VECTORS AND MATRICES* *eeceeecececorseeeerrcerevesssseeseeeneeEe 
hl DIMENSION TRIDIAGONAL MATRICES 


REAL A(N, ¢, NUMCOL), AOLD(N, 4, NUMCOL) 
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=..... DIMENSION STEADY-STATE *t 
REAL SSTI (NUMCOL) 
€ DIMENSION SEMI-INFINITE SOLUTION TIME VECTOR 


REAL TSIS (NUMCOL) 
C.....DIMENSION GEOMETRY VECTOR OF DROPLET HEIGHT 


REAL DEL (NUMCOL) 

DIMENSION SURFACE FLUX VECTOR 

REAL FLUX (NUMNOD) 

DIMENSION LAST SURFACE FLUX VECTOR 

REAL OLDFLX (NUMNOD) 

DIMENSION SURFACE TEMPERATURE VECTOR 

REAL T (NUMNOD) 

DIMENSION THE BEM FORCING FUNCTION TO BE INTEGRATED OVER 

REAL FRCFNC(NUMBET, NUMNOD) 

DIMENSION THE BEM WEIGHT TENSOR USED FOR PRECALCULATED VALUES 

REAL W(NUMNOD, NUMNOD, 10) 

DIMENSION THE VECTOR TO STORE DATA FOR SECOND BEM CALCULATION 

REAL PREDU (NUMNOD) 

DIMENSION THE VECTOR OF SATURATION PRESSURES OF WATER AS FUNCTION OF TEMP. 

REAL PSAT(0:100) 

DIMENSION THE VECTOR OF THE HFG/VFG RATIO FOR WATER AS FUNCTION OF TEMP. 

REAL HVFG(0:100) 

DIMENSION THE VECTOR OF FLAGS FOR INITIALIZATION OF TEMPERATURES OF COLUMNS 

CHARACTER*1 NFTIME (NUMCOL) 

seen DIMENSION THE VECTOR OF TIMES TO OUTPUT THE DATA (.1, .2,..., 3.9, 4., 5., 6.,..., TEND) 
REAL DATTIM(50) 


weerereaeteceneernernetertenenereeeeeeeteraerteseneseeeeeTetentesesteeeetrese 


Sone READ IN TIMES TO OUTPUT DATA (THE NUMBER OF ELEMENTS MUST MATCH THE REAL STATEMENT) 
DATA DATTIM/ 1., 5., 10., 20., 30./ 


suaheteee READ IN PSAT TABLE 
DATA PSAT/.006109, 
+ .006567, .007056, .007577, .008131, .008721, 
+ .009349, .010016, .010724, .011477, .012276, 
+ .013123, .014022, .014974, .015983, .017051, 
+ .018181, .019376, .020640, .021975, .02339, 


+ .02487, .02645, .02810, .02985, -031269. 
+ .03363, .03567, .03782, .04008, .04246, 
+ .04496, .04759, .05034, .05324, .05628, 
+ (059475. " 06281," 066327 7.06999") 7-07304), 
+ .07786, .08208, .08649, .09111, .09593, 
+ 100987 10624, “12275, — (21749) 2349) 
+ 22979, rt3020,,. -14909, 25029, LS 75G, 
+ 7726929, .273931/) “ciSi66 9" 29036, =..2 9940, 
~ ~-20062 el eG0 , a eeR i meee o 84 eee DUS, 
So AeICR Sy -2736, - 2859, - 2986, =P De, 
73296), a9 Si, -3546, 3699, -3858, 
+ .4022, -4192, 4368, .4550, 039. 
- .4934, - S136, 5345, -5560, -5783, 
+ .6014, 6252, - 6498, -6752, -7014, 
+ .7284, . 7564, . 7652, -8149, 8455, 
ane REAP I . 9097, . 9433, .9778, 1.01325/ 


Seo eats READ IN HVFG TABLE 
DATA HVFG/12.1264, 
* 12.9767, 13.8792, 14.8355, 15.8487, 16.9224, 
+ 18.0581, 19.2596, 20.5308, 21.8733, 23.2915, 
+ 24.7897, 26.3694, 28.0369, 29.7935, 31.6445, 
+ 33.5952, 35.6474, 37.8062, 40.0786, 42.4658, 
¢ 44.9764, 47.6111, 50.3778, $3.2823, 56.3274, 
+ 59.5199, 62.8685, 66.3741, 70.0449, 73.8911, 
+ 77.9136, 62.1166, 86.5191, 91.1178, 95.9191, 
+ 100.932, 106.171, 111.638, 117.334, 123.281, 
+ 129.479, 135.931, 142.656, 149.666, 156.964, 
+ 164.551, 172.439, 180.646, 189.182, 198.047, 


102 


"my 


99 


+ 394.431, 410. 
+ 481.498, 500.734, 520.567, 540.948, 562.086, 


+ 583.816, 606. 


+ 1001.597, 


237 373, 247.456, 


207.271, 216.823, 226.778, 


+ 258.877, 270.378, 282.332, 294.671, 307.497, 


320.763, 334.509, 348.725, 363.464, 378.665, 
756, 427.570, 444.927, 462.966, 


207, 629.439, 653.284, 677.668, 
703.116, 729.257, 756.202, 783.711, 812.178, 
841.662, 871.739. 902.695, 934.901, 967.472, 
1035.943, 1071.804, 1108.163, 1146.068, 
1164.769, 1224.430, 1265.227, 1307.048, 1349.996/ 


INPUT USER VARIABLES 
PRINT *, ‘INPUT TEND IN S:’ 
READ (6,*) TEND 
PRINT *, ‘INPUT VO IN MICROLITERS:’ 
READ (6,*) VO 
vo = VO * 1E-09 
PRINT *, ‘INPUT BETAO:’ 
READ (6,*) BETAO 
PRINT *, ‘INPUT THE GEOMETRY MODEL DESIRED, A OR B:’ 
READ (6,*%) GEOMOD 
IF (GEOMOD .EQ. ‘B’) THEN 
PRINT *, ‘INPUT THE INITIAL CONTACT ANGLE THETAO IN DEGREES:’ 


PRINT *, ‘IT MUST BE >= THE MODEL A INITIAL ANGLE: ', CTHETA(BETAQ) * (180. /PI) 
READ (6,*) THETOD 
END IF 


PRINT *, ‘INPUT TSO IN C:’ 

READ (6,*) TSO 

PRINT *, ‘INPUT TCOIL IN C (475-650) :° 

READ (6,*) TCOIL 

PRINT *, ‘INPUT NEGRAT, RATIO USED TO STOP MARCH BACK IN TIME (e.g. 0.01):’ 

READ (6,*%) NEGRAT 

PRINT *, ‘INPUT ROOT NAME OF OUTPUT FILES -- 6 CHARS (e.g. 2run0l):’ 

PRINT *, ‘NAME OF MAIN DROPLET PARAMETERS IS (FILNAM.dat)’ 

PRINT *, ‘NAME OF SOLID SURFACE TEMPERATURES AND FLUXES IS (FILNAM.out)’ 

PRINT *, ‘NAME OF SPATIAL DROPLET SURFACE TEMPERATURES AND FLUXES IS (FILNAM. ins)’ 

PRINT *, ‘NAME OF SPATIAL DROPLET HEIGHT AS FUNCTION OF TIME IS (FILNAM.geo)’ 

PRINT *, ‘NAME OF CUMULATIVE ENERGY TRANSFERS [S -FILNAM.bal)’ 

PRINT *, ‘GIVE OUTPUT FILE NAME INSIDE SINGLE QUOTES IF IT HAS ILLEGAL CHARACTERS’ 

READ (6,°) FILNAM 

PRINT *°, ‘INPUT THE NUMBER OF TIMES TO OUTPUT THE SPATIAL DATA (INTEGER) :’ 

READ (6,°) NDATITM ; 

PRINT ©, ‘INPUT THE TIMES TO OCTPUT THE SPATIAL CATA . 1..2.....3.9,4.,5.,...,TEND): 

30 99 L «© 1, NDATIM 
READ (6,°) DATTIM(L) 

TONT INUE 


CONCATENATE THE REQUIRED OUTPUT FILE NAMES 
FILNAl = FILNAM//’ dat’ 
FILNA2 © FILNAM//' out’ 
FILNAJ e FILNAM//' ins* 
FILNA4 = FILNAM//' geo’ 
FILNAS « FILNAM//' bal’ 


[PEN FILES FOR OUTPUT OF RESULTS 
SPEN (41. FILE «© FILMAL) 
[PEN ‘42, FILE © FILMA?) 
[PEN (43, FILE © FILMA)) 
[PEN ‘44. FILE © FIUNAG) 
CPEN (45. FILE e FILNAS) 


INITIALIZE 
TSFLAG e ' 
TRNBGN « 
EVFLAG e 
THFLAG « 


ooond 
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REMFLG = ‘7’ 


(ag CONVERT TO RADIANS 
THETAO = (PI / 180.) * THETOD 
THETAR = (PI / 180.) * THETRD 


‘ASSUMING MODEL A THE INITIAL THETAA (DEG) = ', CTHETA(BETAO) * (180. / PTI) 


PRINT *, 

Cierarers SET VOLUME TO INITIAL VOLUME 
v = VO 

Corer CALCULATE THE INITIAL RADIUS 


RADO = BETAO * (6. * VO / PI) ** (1. / 3.) / 2. 
RADIUS = RADO 


IF (GEOMOD .EQ. ‘A’) CALL UPGEOA(V, RADO, NUMCOL, DEL) 


IF (GEOMOD .EQ. ‘B’) THEN 
IF (THETAO .LT. CTHETA(BETAO)) THEN 
PRINT *, ‘FOR GEOMETRY MODEL B, THETAO MUST EXCEED THE INITIAL THETA FOR MODEL A!’ 


STOP 
END IF 
CALL UPGEOB(‘1’, V, BETAQ, VO, RADO, THETAO, THETAR, NUMCOL, DEL, RADIUS, THETA) 
END IF 
Cora WRITE THE INITIAL SHAPE TO THE GEOMETRY FILE (*.geo) 


WRITE (44, ‘(13G15.6)‘) TIME, DEL 
PRINT *, ‘SHAPE DISTRIBUTION: ’ 
DO 100 J = 1,NUMCOL 
PRINT *, DEL(J) * 1000. 
100 CONTINUE 


PRINT ts 
PRINT *, ‘RADO:  RADO * T0007 * MM’ 
Soicc ASSUME INITIAL (20C) DENSITY OF WATER IN VOLUME CALCULATIONS 
RHOH20 = 998.2 
Caeramac USE KL AND ALPHA VALUES AT AVERAGE OF TL AND TC FOR SEMI-INFINITE SOLUTION 


PRINT) Sr. ol oUe pla ge las 

TC = SISTC(TSO, TL) 

PRINT *, ‘THE SEMI-INFINITE SOLUTION BOUNDARY TEMPERATURE IS:’, TC 
KLSIS = KL((TL + TC) / 2.) 

ALFSIS = ALPHAL((TL + TC) / 2.) 


Gar ideans SET RADIATIVE PARAMETERS 

EPSF = .2261 

SIGMA = S.67E-08 

ee oFEEF! 

EPSILN = .84 

FRADO = EPSF © SIGMA * (TCOIL + 273 15) °* 4 

FRADS = F * EPSILN © SIGMA © (TCOIL + 273.25) °° 4 

HR = EPSILN * SIGMA © ((TSO + 273.15 - 15.) ¢ 

+ CTAMB (65273515), © ((TSORe 203 ko Lon eu. eee miCAME) ue 62.329) ) ee 2)) 


te Pome aus CALCULATE HCONV, FLUXO (INITIAL CONSTANT FLUX) 
H = HCONV(TSO) 
HO = H + HR 
FLUXO = -SOLVAP(TSO, HO, FRADS) 
PRINT *. ‘FLUXO: ‘, FLOXO 


. INITIALIZE DROPLET TIMES OF SEMI-INFINITE SOLUTION AND NFTIME FLAGS 
-AND FIND THE NUMBER OF TRANSIENT COLUMNS INITIALLY BASED ON THE FOURIER NUMBERS 
DO 110 J = 1, NUMCOL 

TSIS(J) © (DEL(J) °* 2) / (17 © ALFSIS! 


an 


PRINT 6 TSS Go oo ot eee eal 
FO « ALFSIS * DTSHRT / DEL(J) °* 2 
PRINT (ee FON Sarai oop te ist i) 


IF (FO.LT.0.421) NTRANC = J 
NFTIME(J) = ‘0’ 
110 CONTINUE 
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Berge bas INITIALLY MUST USE SHORT TIME INTERVAL 


DT = DTSHRT 
CALL WEIGHT(10, NUMNOD, DT, RADO, W) 


Cereus CALCULATE THE FORCING FUNCTION AT TIME STEP 0.5 
Gos wsiots AND INITIALIZE THE SOLID SURFACE TEMPERATURES 
KS = 1.297 


FLUXOS = KLSIS * (TC - TL) / (PI * ALFSISMS9DY YY i2) ** on5 
DO 120 J = 1,NUMCOL 
T(J) = TC 
FRCFNC(1, J) = -(FLUXOS - FRADO * FRAD(DEL(J), TCOIL) + FLUXO) / KS 
120 CONTINUE 
DO 130 J = NUMCOL + 1 , NUMNOD 
T(J) = TSO 
130 CONTINUE 


Ciara oie CALL BEM1 IN ORDER TO FIND T(J) AFTER ONE TIME STEP 

CALL BEM1(FRCFNC, W, DT, RADO, NUMNOD, NUMBET, NEGRAT, TSO, T, MEMORY, PREDU) 
Career MUST FIND CLOSE VALUES FOR THE INITIAL SSTI(J) VALUES FOR SS COLUMNS 
ase A THE ONE-TERM LINEARIZED TAYLOR SERIES USED CANNOT CORRECT FOR BIG DELTA(SSTI(J)) 


DO 131 J = NTRANC + 1, NUMCOL 
SSTI(J) = TL 
i32 CONTINUE 
SSTI(J) = SSTI(J) + 0.1 


IF (LVFLUX(SSTI(J), PSAT, H, DEL(J), TCOIL, FRADO) - (T(J)-SSTI(3))/DEL(J) .GE. 9.9) GOTO 133 
GOTO 132 
133 CONTINUE 
PRINT *,, “INITIAL SS SSTI(’, J, °):*, SSTI(J) 


131 CONTINUE 


Cetaatere TIME STEP LOOP 
240 CONTINUE 
Ge ocree. IF TIME >= TSHORT THEN WANT TO USE LARGER TIME STEP 
IF (((TIME + 0.0001) .GE. TSHORT) .AND. (TSFLAG .EQ. ‘0’)) THEN 
DT = DTLONG 
TSFLAG = ‘1’ 
CALL WEIGHT(2, NUMNOD, DT, RADO, W) 
CALL RECONF (FRCFNC) 
END IF 


Coren aenecene STEP TIME (TO TIME CORRESPONDING TO END OF THIS ITERATION) 
TIME = TIME + DT 


Sx renctetee RESET CUMULATIVE CHANGE IN VOLUME 

DV = O. 
Cars cieies USE TRIDIAGONAL GAUSSIAN ELIMINATION FOR TRANSIENT COLUMNS 
Corey arars oe8 ONLY IF THE HEAT WAVE HAS REACHED THE LIQUID-VAPOR INTERFACE 

DO 150 J = 1,NTRANC 
a scariest DEL(J) = 0 WAS ALREADY TAKEN CARE OF IN PREDICTOR LOOP 

IF (TIME .GT. TSIS(J)) THEN 
TRNBGN = ‘1’ 


IF (NFTIME(J) .EQ.°0') CALL INTCOL(J., N, TL, T(J), AOLD) 
PRIN ec COG ete me ned Ces) 
NFTIME(J) = ‘1’ 


CALL LINBC(J, N, PSAT, HVFG. H, TCOIL, FRADO. DEL, AOLD, DZ, A) 
CALL GAUSEL(J, N, DEL. DOT. T. A, AOLD, D2) 
SYS onset SAVE TEMPERATURES IN AOLD FOR NEXT ITERATION 


DO 160 I = 1,N 
AOLD(I, 4, J) = A(I, 4. J) 
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n 


1690 CONTINUE 


Ce ecausualiev sre ae PRINT OUT TEMPERATURE COLUMN -- TRANSIENT 
DO 170 I = N,1,-1 
PRINT *, A(I, 4, J) 
170 CONTINUE 
PRINS ee 


CALCULATE SOLID SURFACE FLUXES 
TI = (A(N, 4, J) + A(N - 1, 4, J)) / 2. 
IF (TI .GT. 99.9) THEN 
NFLUX = (KL(100.) * (A(N - 1, 4, J) - A(N, 4, J)) / DZ - H * (100. - TAMB) + 


+ FRADO * (1. - FRAD(DEL(J), TCOIL))) / (HFG(100.) * 18.0152) 
ELSE 
NFLUX = NH20(TI, PSAT, H) 
END IF 


FLUX(J) = KL(T(J)) * (A(l, 4 J) - A(2, 4, J)) / DZ FRADO * FRAD(DEL(J), TCOIL) 


PRINT *, ‘BEM FLUX:’, FLUX(J) 


DVCOL = 18.0152 * DT * NFLUX * PI * RADO ** 2 * (2. * J - 1.) / (NUMCOL ** 2 * RHOH20) 
DV = DV + DVCOL 
PRINT *, ‘TRANSIENT COLUMN DVCOL:’, DVCOL 


ELSE 
Sarat tievaretetete USE THE SEMI-INFINITE SOLUTION FLUX 
PRINT *, ‘HEAT WAVE STILL TRAVELLING IN COLUMN: nia) ed) 
FLUX(J) = KLSIS * (TC - TL) / (PI * ALFSIS * TIME) ** .S - FRADO * FRAD(DEL(J), TCOIL) 
PRINT, (#75 SISFLUX(", J, vi: ', FLUX(J) : 
END IF 
150 CONTINUE 
PRINT *) oe 


Bes SOc UPDATE INTERFACIAL TEMPERATURES FOR STEADY-STATE COLUMNS 
DO 180 J = NTRANC + 1 , NUMCOL 


Sane event IF SOLID IS DRY THEN USE SOLID-VAPOR BOUNDARY CONDITION 
IF (DEL(J) .EQ. 0.) THEN 
FLUX(J) = SOLVAP(T(J), HO, FRADS) 
GOTO 180 
END IF 


IF (SSTI(J) .GE. 100.) SSTI(J) « 99.9 

SSFI «© LVFLUX(SSTI(J), PSAT, H. DEL(J), TCOIL, FRADO) 
SSDFDT = DFDT(SSTI(S). PSAT, HVFG. H) 

SSA =» SSDFDT 

SSB = SSFI - SSTI{J) * SSDFOT 


EXPR = -(SSA * T(J) + SSB) ‘’ ‘SSA * DELIJ) «+ 1.) 
SSTOLD = SSTI(J) 
PRINT ©, “OLD SSTI(°. S01 9" F SSTICI) 


SSTI(J) © DEL(J) *° EXPR + TiJ) 
IF (SSTI‘J) .GT.99.9) THEN 
SSTT(J) © 99.9 
KAVG © KL((T(J) © 100.) / 2.) 
QC @« KAVG © (T(J) - 100.) / DELI(J) 
FLUX(J) = QC - FRADO © FRAD(DEL‘S), TCOIL) 
NFLUX = (QC - H * (100. - TAMB) +¢ FRADO © (1. - FRAD(DEL(J), TCOIL))) / (HFG(100.) * 18.0152) 
PRINT °, ‘SSTI HAS REACHED 100°''''” 
PRINT °*, KAVG. T(J). DELIJ:. OC. FUCK iS) 


ELSE 
PRINT *, ‘SSTI(’'. J, ‘) *. SSTY(J). “OELTSSTI:’. SSTI(J) - SSTOLD 
KAVG «= KL((T(J) + SSTI(J): 2) 


.USE KL VALUE AT AVERAGE TEMPERATURE 

FLOX(J) = -KAVG * EXPR - FRADO * FRAD(DEL(J). TCOIL) 
NFLUX © NH20(SSTI(J), PSAT. H) 

END IF 

PRINT °, ‘'SSFLUX:’, FLUX(J) 


106 


DVCOL = 18.0152 = DT * NFLUX * PI * RADO ** 2 * (2. © J - 1.) 4 ‘NUMCOL ** 2 


DV = DV + DVCOL 
PRINT *, ‘STEADY-STATE COLUMN DVCOL:’, DVCOL 


CALL UPFFNC(‘1’, FLOX, FLOXO, NUMBET, NUMNOD, FRCFNC) 
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RHCH20) 


180 CONTINUE 
PRINT *, ‘THE LAST TRANSIENT COLUMN IS COLUMN: ', NTRANC 
ao ee CHECK fLUX RATIO OF OUTERMOST TRANSIENT COLUMN TO SEE IF IT SHOULD BE SS 
IF ((NTRANC .GT. 0) .AND. (TRNBGN .EQ. ‘1‘)) THEN 
J = NTRANC 
FLUXR auta(Nee4ee Ties AUNe= dser4y-15) 0) / (A(2, 4, 3) - All, 4,090 
PRINT *, ‘FLUXR(’, J, ‘):', FLUXR 
IF (FLUXR .GT. PERC / 100.) THEN 
SSTI(J) = (A(N, 4, J) + A(N- 1, 4, J)) / 2. 
Celene TRANSFER TI 
NTRANC = NTRANC - 1 
IF (NTRANC .EQ. 0) SSTIME = TIME 
END IF 
END IF 
cohras semes UPDATE THE VOLUME OF THE DROPLET 
VeV- DV 
PRINT *, ‘VOLUME:’, V * 1E+09, ’ MICROLITERS’ 
IF ((V .LT. 0.) .AND. (EVFLAG .EQ. ‘0’)) THEN 
PRINT *, ‘THE COLUMN IS COMPLETELY EVAPORATED. ' 
DO 190 J = 1,NUMCOL 
DEL(J) = 0. 
190 CONTINUE 
EVFLAG = ‘1’ 
EVTIME = TIME 
END IF 
IF (EVFLAG .EQ. ‘0’) THEN 
IF (GEOMOD .EQ. ‘B’) THEN 
CALL UPGEOB(’0’, V, BETAO, VO, RADO, THETAO, THETAR, NUMCOL, DEL, RADIUS, THETA) 
IF ((THETA .LE. THETAR) .AND. (THFLAG .EQ. ‘0’)) THEN 
THTIME = TIME 
THFLAG = ‘1’ 
END IF 
END IF 
Cres Cae CALCULATE THE CURRENT BETA (BASED ON CURRENT V AND RADIUS) 
BETA = 2. * RADIUS / (6. * V / PI) °* (1./3.) 
IF (GEOMOD .EQ. ‘A’) THEN 
CALL UPGEOA(V, RADO, NUMCOL, DEL) 
THETA = CTHETA(BETA) 
END IF 
PRINT *, ‘SHAPE DISTRIBUTION: ’ 
DO 200 J = 1,NUMCOL 
PRINT ©, DEL(J) © 1000. 
200 CONTINUE 
PRINT ©, ‘ ° 
END IF 
ere. CALCULATE SURFACE FLUXES FOR R > RADO 
DO 220 J = NUMCOL + 1 . NUMNOD 
FLUX(J) = SOLVAP(T(J), HO, FRADS) 
220 CONTINUE 
Date 6. SAVE ALL SURFACE FLUXES IN OLDFLX() 
DO 230 J = 1,NUMNOD 
OLDFLX(J) © FLUX(J) 
230 CONTINUE 


201 


207 


204 


206 


202 


WRITE DATA TO OCUTPUT FILES weereerverereeereececrevereerestvereverrentereewereTereT 


PRINTS 2a 
PRINT *, ’ J FLUX (J) T(J) TI SSTI’ 
DO 201 J = 1, NUMCOL 
TICONV = (AOLD(N, 4, J) + AOLD(N - 1, 4, J)) / 2. 
PRINT *, J, FLUX(J), T(J), TICONV, SSTI(J) 
CONTINUE 
DO 206 L = 1, NDATIM 
IF (INT(10.*DATTIM(L)) .EQ.INT(10.*TIME+0.01)) THEN 
DO 207 J = 1, NUMNOD 
IF (J .LE. 48) THEN 
R = RADO * (J - .5) / 12. 
ELSEIF (J .LE. 66) THEN 
R = RADO * (4. + ((J - 48) - .S) / 6.) 
ELSEIF (J .LE. 75) THEN 
R = RADO * (7. + ((J - 66) - .S) / 3.) 
ELSE 
R = RADO © (10. + ((J = 75) = <25)) 
END IF 
WRITE (42, ‘(6G14.6)’) TIME, R, R / RADO, T(J), FLUX(J), FRCFNC(1,J) 
CONTINUE 
DO 204 J = 1, NUMCOL 
IF (TIME.LE.TSIS(J)) THEN ~° 
TICONV = TL 
ELSEIF (J.LE.NTRANC) THEN 
TICONV = (AOLD(N, 4, J) + AOLD(N - 1, 4, J)) / 2. 
ELSE 
TICONV = SSTI(J) 
END IF 
WRITE (43, ‘(G14.6, I4, 4G14.6)‘) TIME, J, FLUX(J), FRCFNC(1,J), T(J), TICONV 
CONTINUE 
END IF 
CONTINUE 
WRITE (44, ’(13G15.6)’) TIME, DEL 
PRINT *, ‘FLOXR: ‘, FLUXR 
PRINT *, “VOLUME: ‘, V * 1E+09, °TIME: ", TIME, “ S* 
PRINT *, ‘DV: ’, DV ® 1E+09 
PRINT *, ‘MEMORY: ‘, MEMORY 
PRINT *, ‘BETA: ‘, BETA, ‘THETA: ‘, THETA * (180. / PI) 
APPROXIMATE CHECK OF OVERALL ENERGY BALANCE 
TICAVG = 0. 
TAVG = 0. 
FAVG = 0. 
DO 202 J = 1, NUMCOL 
TAVG = TAVG + T(J) © (2. * J - 1.) / NUMCOL ** 2 
FAVG = FAVG + FLUX(J) * (2. * J - 1.) / NUMCOL ** 2 
CONDUC = CONDUC + DT * FLOX(J) * PI * RADO ** 2 * (2. * J - 1.) / NUMCOL ** 2 
RADIAT = RADIAT + DT * FRADO * (1. - FRAD(DEL(J), TCOIL)) * PI * RADO ** 2 * (2. 
IF (TIME.LE.TSIS(J)) THEN 
TICONV = TL . 
ELSEIF (J.LE.NTRANC) THEN 
TICONV = (AOLD(N, 4, J) + AOLD(N - 1, 4, J)) / 2. 
ELSE 
TICONV = SSTI(J) 
END IF 
TICAVG = TICAVG + TICONV * (2. * J - 1.) / NUMCOL ** 2 
CONVEC = CONVEC + DT © H © (TICONV - TAMB) * PI * RADO ** 2 * (2. * J - 1.) 
CONTINUE 
LATENT = LATENT + DV * 998.2 * HFG(TICAVG) 
MAXSEN = V * 998.2 * (440. - 104.89) * 1000. 
SENSIB = SENSIB + DV * 998.2 * (440. - 104.89) * 1000. 
PRINT *, ‘CONDUC: ‘, CONDUC, ‘RADIAT: ’, RADIAT 
PRINT *, ‘CONVEC: ’, CONVEC, ‘MAXSEN: ', MAXSEN, ‘SENSIB: ‘, SENSIB 
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J 


1.) / NUMCOL ** 2 


/ NUMCOL ** 2 


240 


PRINT ©, ‘TICAVG: ', TICAVG, ‘TAVG: ‘, TAVG, LATENT: ‘, LATENT 
PRINT *, ’BALANCE: ‘, CONDUC + RADIAT - CONVEC - LATENT 
WRITE (41, '(10G14.6)’) TIME, FLUXR, V * 1E9, DV * 1E9 / DT, MEMORY, 
BETA, THETA * (180. / PI), TAVG, TICAVG, FAVG 
WRITE (45, ‘(7G14.6)’) TIME, CONDUC, RADIAT, CONVEC, LATENT, SENSIB, MAXSEN 


rrr rrrrrTrT TTT Titi titties 


CALL BEM1(FRCFNC, W, DT, RADO, NUMNOD, NUMBET, NEGRAT, TSO, T, MEMORY, PREDU) 


$44444444+4++4++4NOW PREDICT THE NEXT FLUX VECTOR+e++++esrerersosssoses 
PRINT *, ‘NOW GOING INTO PREDICTION LOOP!!!!!’ 
DO 240 J = 1, NTRANC 
IF COLUMN EVAPORATED BEFORE STEADY-STATE THEN TRANSFER TO SS 
IF (DEL(J) .£Q. 0.) THEN 
NTRANC = J - 1 
GOTO 250 
END IF 
IF (TIME .GT. TSIS(J)) THEN 
CALL LINBC(J, N, PSAT, HVFG, H. TCOIL, FRADO, DEL, AOLD, DZ, A) 


CALL GAUSEL(J, N, DEL, DT, T. A. AOLD, DZ) 


eee Ae DO NOT SAVE TEMPERATURES IN AOLD FOR NEXT ITERATION 


sn a PRINT OUT TEMPERATURE COLUMN -- TRANSIENT 


PRINT *, ’THESE TEMPERATURES ONLY USED TO PREDICT FLUX()’ 
DO 260 I = N,1,-1 

PRINT *, A(Z, 4, J) 
CONTINUE 


oterate ts CALCULATE PREDICTED SOLID SURFACE FLUXES FOR R < RADO 
FLOX(J) = KL(T(J)) © (All, 4, J) - Al2, 4, J)) / DZ - FRADO © FRAD(DEL(J), 


PRINT *, ‘BEM PREDICTED FLUX:’, FLUX(J) 
ELSE 


Meriess USE THE SEMI-INFINITE SOLUTION FLUX 


TCOIL) 


FLUX(J) = KLSIS * (TC - TL: © ‘PI * ALFSIS ° TIME) ** -.5 - FRADO * FRAD(DEL(J), 


PRINT *, ‘SISFLOX('. J. ‘} ‘, FLOX.J) 
END IF 
CONTINUE 
PRINS 


CONT INUE 


.. CONTINUE PREDICTION OF FLUXES WOW :M THE STEACY-STATE MODE 


DO 270 J = NTRANC + 1 . NUMCCL 


...IF SOLID IS DRY THEW USE SOLID-VAPOR BOCNDARY CONTITION 
IF (DEL(J) .E—Q. 0.) THEW 
FLUX(J) © SOLVAP(T(J). HO. FRADS) 
GOTO 270 
EXD IF 


IF 'SSTI(J) .GE. 100 ) SSTI‘J' o 99 9 

SSFI e LVFLOX(SSTI(J). PSAT. MH OEL SL TCSlL. ~FRADO) 

SSDFOT «© OFDT(SSTI iJ). PSAT. MVPG A 

SSA « SSDFDT 

SSB «© SSFI - SSTI(J) ° SSOPUOT 

EXPR © - (SSA * TiJ) + SSB) S$A ° DELS +e i 

-0DO0 MOT UPDATE THE STEADY-STATE 2: 5Ci5-VAPCR INTERFACIAL TEMPERATURE 
SSTEMP © DEL(J) * EXPR - TiJ 

IP (SSTEMP GT.99 9: THEW 


RAVG © KL((T'J) « 100) 2 

QC e KAVG * (TiJ) - 100 cE. ZL 

FLOXiJ) « QC - FRADO © FRAD CEL L TCCIL 
ELSE 

BAVG e KL((T(J) © S8TI(J)) » 32 
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TCOIL) 


S OLD KAVG VYALTE GOOD ENOUGH 


FLUX(J) = -KAVG * =XPR - FRADO * FRAD(DEL(J). TCOIL) 


END IF 
PRINT *, ‘PREDICTED SSFLUX(’, J, ‘):'. FLUX(J) 


270 CONTINUE 


(abe Sas eea CALCULATE SURFACE FLUXES FOR R > RADO 
DO 280 J = NUMCOL + 1 , NUMNOD 
FLUX(J) = SOLVAP(T(J), HO, FRADS) 
280 CONTINUE 


Cc eee we PEESS ESSE SSE SEES HHH S HSH SHH HSH HS SSH HHH FHS HHH SHEP HHH HP SP PHP HEHEHE HE+OHH HH 


feoscnAF AVERAGE FLUX AND OLDFLX TO GET FLUX AT CENTER OF CURRENT TIME INTERVAL 
DO 290 J = 1,NUMNOD 
FLUX(J) = (FLUX(J) + OLDFLX(J)) / 2. 
290 CONTINUE 


CALL UPFFNC('0’, FLUX, FLUX0, NUMBET, NUMNOD, FRCFNC) 
CALL BEM2(FRCFNC, W, NUMBET, DT, RADO, NUMNOD, TSO, T, PREDU) 
IF ((TIME + 0.0001) .GE. TEND) GOTO 9999 
GOTO 140 
9999 CONTINUE 
Caer PRINT OUT RESULTS 
PRINT *, ‘TIME UNTIL THE DROPLET EVAPORATED = ', EVTIME, ‘S’ 
PRINT *, ‘TIME UNTIL QUASI-STEADY STATE WAS REACHED = ‘, SSTIME, ‘S’ 


IF (GEOMOD .EQ. ‘B’) PRINT *, ‘TIME UNTIL THETA REACHED THETAR WAS = ', THTIME, ‘S’ 


STOP 
END 


REAL FUNCTION ALPHAL (T) 


Giaenens ALPHAL THERMAL DIFFUSIVITY OF WATER (M**2/S) 

(SES 95.5 USE FIFTH ORDER POLYNOMIAL TO APPROXIMATE ALPHA FOR WATER AS FUNCTION OF TEMPERATURE 
REAL T 
ALPHAL? = 0000002) 29123514 20082737 = Tl= 3) 3096E-05 * “h ** 25+ 3. 96S6E-07 2 T9253 
+ = 3. SESE-O09)"* | 7 eRe ei issy Rink Lets Tae? 15)) 
RETURN 
END 


SUBROUTINE BEM1 (FRCFNC, W, DT. RADO, NUMNOD. NUMBET. NEGRAT, TSO. T, MEMORY, PREDU) 


iSeacas PERFORM BOUNDARY ELEMENT INTEGRATION USING WEIGHT TENSOR AND FORCING FUNCTION 
Cosinm Ue T + *tS6 
Caen UMAX =» MAX ADDITION TO U AT TIME TO 


INTEGER NUMNOD, NUMBET 

REAL FRCFNC(NUMBET, NUMNOD), W(NUMNOD. NUMNOD, 10). T(NUMNOD). PREDU (NUMNOD) 
REAL NEGRAT, TSO, DT, RADO, MEMORY 

REAL U(78), UMAX, UTERMK 

INTEGER I, J, K 

CHARACTER*®1 EXFLAG (786) 


Caeipar U: ACCUMULATING TEMPERATURE DEPRESSION: UMAX MAGNITUDE OF LARGEST UTERMK; 
Gai srarys UTERMK: TERM IN BEM INTEGRATION; 
Caererata EXFLAG: FLAG -- ‘1’ » NODE READY TO IT RECOLLECTION TIME LOOP 

REAL WGHT 


UMAX = 1.E-10 
DO 400 J = 1,NUMNOD 
Utd) = O. 
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PREDU(J) = 0. 
EXFLAG(J) = ‘0’ 
400 CONTINUE 
DO 410 K = 1,NUMBET 
DO 420 J = 1,NUMNOD 
UTERMK = 0. 
DO 430 I = 1,NUMNOD 
UTERMK = UTERMK + WGHT(J, I, K, DT, RADO, W) * FRCFNC(K, I) 
430 CONTINUE 
U(J) = U(J) + UTERMK 
IF (K.GT.1) PREDU(J) = PREDU(J) + UTERMK 
IF (ABS(UTERMK) .GT. UMAX) UMAX = ABS (UTERMK) 
IF (ABS(UTERMK) / UMAX .LE. NEGRAT) EXFLAG(J) = ‘1’ 
420 CONTINUE 
DO 440 J = 1,NUMNOD 
IF (EXFLAG(J) .EQ. ‘'0’) GOTO 410 
440 CONTINUE 
GOTO 450 
410 CONTINUE 
PRINT *, ‘NOT ENOUGH TIME STEPS IN THE INTEGRATION WERE USED!!!’ 
STOP 
450 CONTINUE 
DO 460 J = 1,NUMNOD 
T(J) = U(J) + TSO 
460 CONTINUE 
MEMORY = K * DT 
PRINT *, ‘THE NUMBER OF BEM TIME STEPS USED WAS: eae 4 
RETURN 
END 


SUBROUTINE BEM2 (FRCFNC, W, NUMBET, DT, RADO, NUMNOD, TSO, T, PREDU) 
ai niee PERFORM BOUNDARY ELEMENT INTEGRATION USING WEIGHT TENSOR AND FORCING FUNCTION 
crits CALCULATIONS ARE SHORTENED BY USING PREDU FROM BEM1 
= highs U = T - TSO 

INTEGER NUMNOD, NUMBET 

REAL FRCFNC(NUMBET, NUMNOD), W(NUMNOD, NUMNOD, 10), T(NUMNOD), PREDU(NUMNOD) 

REAL TSO, DT, RADO 

REAL UTERMK 

INTEGER I, J, K : 
ee iat UTERMK: TERM IN BEM INTEGRATION 

REAL WGHT 
25 Oe ONLY NEED TO RECALCULATE THE FIRST TIME INTERVAL BACK IN TIME 

Kel 

DO 470 J = 1,NUMNOD 

UTERMK e 0. 
DO 480 I = 1,NUMNOD 
UTERMK = UTERMK + WGHT(J, I. K, OT, RADO, W) © FRCFNC(K, I) 
480 CONTINUE 
T(J) = PREDU(J) + UTERMK + TSO 
470 CONTINUE 
RETURN 
END 


REAL FUNCTION CTHETA (BETA) 
CALCULATE THE CONTACT ANGLE GIVEN BETA ASSUMING GEOMETRIC MODEL A (SEGMENT OF SPHERE) 
REAL BETA, GAM 
Sear GAM: THICKNESS OF DROPLET AT APEX / RADO 
GAM = (4./BETA ** 3 + (1. + 16./BETA °°6) °*.5) ** (1./3.) 


- - (+4. /BETA *°3 ¢ (1. ¢ 16./BETA **6) °¢.S) ** (1./3.) 
CTHETA © ATAN(((1./GAM + GAM) °°2 /4@. - 1.) ** -.5) 

RETURN 

END 
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REAL FUNCTION CPA {T) 


_CPA CONSTANT PRESSURE SPECIFIC HEAT FOR AIR (J/KG-K) 
SSE SECOND ORDER POLYNOMIAL TO APPROXIMATE CP FOR AIR AS FUNCTION OF TEMPERATURE 


REAL T 

CPA = 1003.4 + .031646 * T + 3.4286E-04 * T ** 2 
RETURN 

END 


REAL FUNCTION DFDT (TI, PSAT, HVFG, H) 

CALCULATE DERIVATIVE OF YAPOR-LIQUID FLUX WRT TEMPERATURE 

REAL TI, PSAT(0:100), HVFG(0:100), H 

REAL TAMB, TF. XA, PATM, LEZ3, COEFF, XTERM 

TAME: AMBIENT TEMP; TF: BOUNDARY LAYER FILM TEMP; XA: AMEILENT MOLE FRACTION OF WATER; 


PATM- ATMOSPHERIC PRESSURE: LE23- LEWIS MUMBER * (2/3); 


_COEFF, XTERM: ALGEBRAIC EXPRESSIONS 


S 


PATM = 101325. 

LE23 = .494 

COEFF = .621924 * H * HFG(TI) / (CPA(TF) * LE23) 

XTERM = (1. - XI(TI, PSAT)) ** 2 

DFDT = (COEFF * (1 - XA) * HV(TI, HVFG) / (PATM * XTERM * (TI + 273.15)) + H) / KL(TI) 
CAH USE DFDT = H / KL(TI) AS A CHECK FOR CASE OF CONVECTION ONLY 


RETURN 
END 


REAL FUNCTION FRAD (D, T) 

CALCULATE FRACTION OF RADIATIVE FLUX REACHING DEPTH D IN DROPLET 
FUNCTION IS A 2D CURVE FIT OF RADIATION STUDY RESULTS 

REAL D, T 

REAL Z. MK, B 


-D: DEPTH BELOW LIQUID-VAPOR INTERFACE (mM) 


2: DEPTH (M4) 


I: TEMPERATURE OF RADIATIVE HEATER COILS {C) 


FRAD: MORMALIZED RADIATIVE FLUX (FLUX/[(EPS*F) *SIGMA*® (T+273.15) **4)) 
Ze<«D* 1000. 
IF (Z .LE. .04) THEN 

Ms 322.37 - .$6492 * T + 3.2645E-04 * FT °° 2 

Bs 0. 

ELSEIF ‘2 -LE. .1) THEN 

Me 1042.3 - 2.7002 * T + 0018666 * T #9 

BS © -61.573 + .29781 © F - 3.6576E-06 © T 9° 2 ©: 1.7732E-07 © T 29 3 
ELSEIF (Z LE. .2) THES 

MS © 3863.7 - 16-516 © FT e .626657 © FT oe 2 - 1.2372E-05 * FT #¢ 3 

Boe -273. + 1.2567 © F - .0039633 ° T #* 2 + 1.00968-06 *° T *° 3 
ELSEIF ‘2 LE. .7) THEW 

M © 4122.1 - 18-206 * FT + .027637 * FT °* 2 - 1. 4258B-05 * T 2? 3 

Be -309.61 + 1.95276 © F - 06634716 ° F °° 2 + 1.33038-06 ° T ** 3 
ELSE 


6S « $577.9 + 26:76 8 FT 0° /OSTI03 OF Kee 2 - 3, 9699E-OS.° T.9*.3 
Bo -1336. + 6.14869 © FT - -90095698 ° 7 9° 2 + $.02268-06 * T *9° 3 
EOD IF 
FRAD © 1. / (© Ze Be 1.) 
DETTRN 
L > $9) 


SUBROUTINE GAUSEL (J. MW, DEL, OT. T. A. AOLD, DZ) 

USE TRIDIAGOMAL GAUSSIAN ELIMIMATION TO UPDATE TEMPERATURES IN DROPLET COLUMN 
PAZSKETER (HUMCOL «= 12) 

INTEGER J, W 

BEAL DEL (WUMCOL). OT, TIMUMCOL). A‘l2. 4, WOMCOL), AOLD(i2, 4, MUMCOL), DZ 
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L Ti. ALPHA. SAMMA 


TZ2: CRANK-NICHOLSON EXPRESSION: ALPHA: THERMAL DIFFUSIVITY CF LiQurin: 


GAMMA: CRANK-NICHOLSON CONSTANT 
REAL ALPHAL 
Ati, 1. J) = 0. 
A(l, 2, J) = .5 
ati. 3, J) = .5 
A(i, 4, J) = Tg) 
Tz2 = UT f/f t2~* DE **' 2) 
THE C-N ALPHA IS A MEASURE OF ACCURACY OF 
GIVE A CONSERVATIVE (UPPER LIMIT) VALUE 
PRINT *, ‘C-N ALPHA(’, J. "): ag a; uk 
bo S00. 2:.<°2°,. 8 ="t 
ALPRA = ALPHAL(AOLD(I, 4, J)) 
GAMMA = ALPHA * TZ2 
A(I, 1, J) = -GAMMA 
A(I, 2, J) = 1. + 2. * GAMMA 
A(I, 3, J) = -GAMMA 


THE CRANK-NICHOLSON METHOD 


7E-07 © Tz2 


A(I, 4, J) = AOLD(I, 4, J) + GAMMA © (AOLD(I - 1, 4. J) - 2. * AOLD(I, 


CONTINUE 


eererereerereeeMEAT OF GAUSSIAN ELIMINATION* feeererecererrecee 


DO $i0 I = 2,N 


Atl. 2, di) @ Alt, 4; 2)) = ANZ, 2. 9d) 7 AG = 2. 2. J} * AZ = 2, 3. JI 
Ati, 4. 3) = A(I, & J) - Atl, % J) / AlT = 1.. 2. F * AlT =... 4. J) 


CONT. SUE 
A(N, «, J) = A(N, 4, J) / A(N, 2, J) 
bo $20 Te=eN-2,i1,. =2 
A(I, 4. J) = (A(I, 4, J) - AlI, 3, J) * 
CONTINUE 


eeecerrerrererrrrererrererterererrrereeernrrerereereerrererereereres 


RETURN 
END 


REAL FUNCTION KCONV (T) 


A(I + 1,.4, J)) / A(t. 2. J) 


HCONV EXPERIMENTAL CONVECTIVE HEAT TRANSFER COEFFICIENT (WN/M**2-K) 
USE THIRD ORDER POLYNOMIAL FOR EXPERIMENTAL HCONV AS FUNCTION OF TSO 


REAL T 


HCONV © -42.348 + 1.3663 * T - .011498 © T 


RETURN 
END 


REAL FUNCTION HFG (TI) 

CALCULATE THE LATENT HEAT OF WATER (J/KG) 
REAL TI 

REAL T, Al. AZ, A3, Ad. AS, AG. LAMBDA 
Te (647.27 - (TI * 273.35)) # 647.27 
Al = .72241 

A2 = §.33402 

AJ «= 8.97347 

A4@ = -11.93143 

AS » -3.31206 

A6 = 1.63257 


EAGGDA = Al? Tf Stl. 73.) © ad SY ee 29 2 


HFG = LAMBDA * 2501000. 
RETURN 
END 


REAL FUNCTION HV (TI. HVFG) 


ee 2 + 3.1956E-05 © T ef 3 


AS A FUNCTION OF TEMPERATURE 


AD Mel 208 2 AE © 7 2 AS? T 


FUNCTION TO GIVE HFG/VFG GIVEN TEMPERATURE 
HV: RATIO OF LATENT HEAT TO SPECIFIC VOLUME CHANGE OF VAPORIZATION 


REAL TI, HVFG(0:100) 
INTEGER TSAT 
REAL TVAL, HVA, HVB 
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J) + AOLD(I + 2, 


ee2 2 AG * T ee3 


4. 


J)) 


TSAT = -1 

TVAL = TI 

IF (TVAL .GT. 100.) TVAL = 190. 
550 CONTINUE 

TSAT = TSAT + 1 

IF (TSAT .LT. TVAL - 1.) GOTO S50 
HVA = HVFG(TSAT) * 1000. 
HVB = HVFG(TSAT + 1) * 1000. 
LINEARLY INTERPOLATE HV 
HV = HVA + (HVB - HVA) * (TI - TSAT) 
RETURN 
END 


SUBROUTINE INTCOL (J, N, TL, TC, AOLD) 
Setar CALCULATE TSIS(J), THE TIME UNTIL HEAT WAVE IN COLUMN J REACHES SURFACE AND 


INITIALIZE TRIDIAGONAL MATRIX USING SIS SOLUTION TEMPERATURE PROFILE 
PARAMETER (NUMCOL = 12) 
INTEGER J, N 
REAL TL, TC, AOLD(12, 4, NUMCOL) 
REAL ERF 
eS ERF: ERROR FUNCTION (FROM LINKED LIBRARY) 
DO. 600 } = 2, N- 12 

AOLD(I, 4, J) = (TL - TC) * ERF((I - 1.5) * 17. **.5 / (2. * (N - 2.))) + TC 
600 CONTINUE 
ADJUST Tl AND TN IN ORDER TO GET CORRECT BOUNDARY TEMPERATURES 
AOLD (2 04g) n2s 2) TC -eAOLD (2) 41110) 
AOLD(N, 4, J) = 2. * TL - AOLD(IN - 1, 4, J) 
RETURN 
END 


REAL FUNCTION KL (T) 
<esteheee KL THERMAL CONDUCTIVITY OF WATER (W/M-K) 
USE FIFTH ORDER POLYNOMIAL TO APPROXIMATE K FOR WATER AS FUNCTION OF TEMPERATURE 


REAL T 

Ki = .56971 + .001754 * T - 4.0332E-06 * T **2 = 4.3732E-08 * T **3 
+ * 2.202E-10 * T °84 = 2.9455E-13 * T ¢*5 
RETURN 

END 


SUBROUTINE LINBC (J, N, PSAT, HVFG, H, TCOIL, FRADO, DEL, AOLD, DZ, A) 
LINEARIZE THE VAPOR BOUNDARY CONDITION USING FIRST TERM OF TAYLOR SERIES 
PARAMETER (NUMCOL = 12) 
INTEGER J, N 
REAL PSAT(0:100), HVFG(0:100), H, TCOIL, FRADO, DEL(NUMCOL) 
REAL AOLD(12, 4, NUMCOL), DZ, A(12, 4, NUMCOL) 
REAL TIO, FIO, DFDTO, AA, BB 
eer cy 65 TIO, FIO, DFDTO, AA, BB: PARAMETERS FOR TAYLOR LINEARIZATION OF LIQUID-VAPOR BC 
REAL DFDT, LVFLUX 
DZ = DEL(J) / (N - 2.) 


PRINT *, ‘DZ(’, J, ‘):’. DZ * 1000., ‘MM’ 
TIO = (AOLD(N, 4, J) + AOLD(N - 1, 4, J)) / 2. 
PRINT ©, O° TEOs ) TIO, “TUG WS, © ee CAOLD( 145 ot) 2 AOLD(2, 4, J)) / 2: 


IF (TIO .GT. 99.9) TIO © 99.9 
FIO = LVFLUX(TIO, PSAT, H, DEL(J), TCOIL, FRADO) 
DFDTO = DFDT(TIO, PSAT, HVFG, H) 
AA © DZ * DFDTO 

BB = DZ * (FIO - TIO * DFDTO) 
A(N, 1, J) #© 1. - AA / 2. 

A(N, 2, J) © -1. - AA / 2. 

A(N, 3, J) = 0. 

AIN, 4, J) = BB 

RETURN 

END 


114 


non 


a) 


REAL FUNCTION LVFLUX (TI, PSAT, H, D, TCOIL, FRADO) 
Sarto FUNCTION TO CALCULATE FLUX AT LIQUID-VAPOR INTERFACE 

REAL TI, PSAT(0:100), H, D, TCOIL, FRADO 

REAL TAMB, TF, XA, LE23, COEFF, XTERM, CONV, RAD 
Senet TAMB: AMBIENT TEMP; TF: BOUNDARY LAYER FILM TEMP; XA: AMBIENT MOLE FRACTION OF WATER; 
Sis spies LE23: LEWIS NUMBER * (2/3); COEFF, XTERM: ALGEBRAIC EXPRESSIONS; 
wis eee CONV: CONVECTIVE HEAT TRANSFER; RAD: RADIATIVE HEAT TRANSFER 

REAL HFG, CPA, XI, FRAD, KL 

TAMB = 25S. 

TF = (TI + TAMB) / 2. 

XA = O. 

LE23 = .894 

COEFF = .621924 * H * HFG(TI) / (CPA(TF) * LE23) 

XTERM = (XI(TI, PSAT) - XA) / (1 - XI(TI, PSAT)) 

CONV = H * (TI - TAMB) 

RAD = FRADO * (1. - FRAD(D, TCOIL)) 

LVFLUX = (COEFF * XTERM + CONV - RAD) / KL(TI) 

RETURN 

END 


REAL FUNCTION NH20 (TI, PSAT, H) 
Saas CALCULATE MOLAR FLUX AT LIQUID-VAPOR INTERFACE FOR TI < 100 
REAL TI, PSAT(0:100), H 
REAL TAMB, TF, XA. LE23, MA, COEFF. XTERM 
oes TAMB: AMBIENT TEMP; TF: BOUNDARY LAYER FILM TEMP; XA: AMBIENT MOLE FRACTION OF WATER; 
Pa ae LE23: LEWIS NUMBER ~ (2/3): MA: MOLECULAR MASS OF AIR 
aes COEFF, XTERM: ALGEBRAIC EXPRESSIONS 
REAL CPA, XI 
IF (TI .GE. 100.) THEN 
PRINT *, ‘NH20 CALCULATION TRIED FOR TI >= 100. !’ 
STOP 
END IF 
TAMB = 25. 
TF = (TI + TAMB) / 2. 
XA «= 0. 
LE23 = .894 
MA = 28.9669 
COEFF «= H / (MA * CPA(TF! ° LE2) 


XTERM © (XI(TZ. PSAT) - XA) 4 ‘3 - KI(T2. PSAT)! 
NH20 # COEFF ° XTERM 

RETURN 

END 


SUBROUTINE RECONF | FRCFWC) 
FIX THE FORCING FUNCTION FOR THE SWITCH TO STLOMG 
WARNING: GOOD ONLY FOR OTSHORT © © :. OTLOMG @ i 3. AND TSHORT © 4. !:!!!!! 
REAL FRCFNC(100. 78) 
INTEGER J. K 
D0 610 Je 1.78 
dO 620 Ke l.4@ 
FRCFNC(K, J) e (FRCPIC(30°R-$ 5 + FRCOPWC.19°K-4, J): 2 
620 CONTINUE 
30 630 K e $,100 
FRCFNC(K. J) © 0 
630 CONT INUE 
€19 CONTINUE 
RETURN 
END 


REAL FUNCTION SISTC (TSO TL 
CALCULATE THE SEMI-INFINITE SOLCTIOM COMTACT TEMPERATURE OF DROPLET AND SOLID 
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ASSUMING INITIAL TEMPERATURE OF SOLID IS CONSTANT (IT IS ACTVALLY =INEAR) 


REAL TSO, TL 

REAL RHOS, CPS, KS, RHOL, CPL, RADS, RADL, NUM, DEN 

RHOS: DENSITY OF SOLID; CPS: SPECIFIC HEAT OF SOLID; KS: CONDUCTIVITY OF SOLID; 
RHOL: DENSITY OF LIQUID; CPL: SPECIFIC HEAT OF LIQUID; 

RADS, RADL, NUM, DEN: ALGEBRAIC EXPRESSIONS 


REAL KL 

RHOS = 2520. 

CPS = 888.9 

KS = 1.297 

RHOL = 998.2 

CPL = 4179. 

RADS = (RHOS * CPS * KS) ** .5 
RADL = (RHOL * CPL * KL(TL)) ** .5 
NUM = TSO * RADS + TL * RADL 
DEN = RADS + RADL 

SISTC = NUM / DEN 

RETURN 

END 


REAL FUNCTION SOLVAP (T, HO, FRADS) 

CALCULATE THE FLUX FOR DRY SURFACE AT TEMPERATURE T 
SOLVAP IS POSITIVE OUT OF THE SURFACE, SOLVAP < 0 
REAL T, HO, FRADS 

REAL TAMB 

TAMB: AMBIENT TEMPERATURE 

TAMB = 25. 

SOLVAP = HO * (T - TAMB) - FRADS 

RETURN 

END 


SUBROUTINE UPFFNC (RLLDWN, FLUX, FLUXO, NUMBET, NUMNOD, FRCFNC) 
UPDATE THE ARRAY OF TEMPERATURE GRADIENT FORCING FUNCTIONS 
THROW OUT THE OLDEST IF RLLDWN = ‘1’ 
CHARACTER*1 RLLDWN 
INTEGER NUMBET, NUMNOD 
REAL FLUX(78), FLUXO, FRCFNC(100, 78) 
REAL KS 
INTEGER I, J 
KS: CONDUCTIVITY OF SOLID 
KS = 1.297 
IF (RLLDWN .EQ. ‘1°’) THEN 

DO 650 I = NUMBET,2,-1 

DO 660 J = 1,NUMNOD 
FRCFNC(I, J) =» FRCFNC(I - 1. J) 
CONT INUE 

CONTINUE 
END IF 
DO 670 J = 1,NUMNOD 

FRCFNC(1, J) = -(FLOX(J) + FLOXO) / KS 
CONTINUE 
RETURN 
END 


SUBROUTINE UPGEOA (V, RADO, NUMCOL. DEL) 

UPDATE THE HEIGHT OF THE DROPLET USING MODEL ‘A’ 
PARAMETER (PI © 3.14159265358979) 

REAL V, RADO, DEL(12), BETA 

INTEGER NUMCOL 


REAL GAM, R 
GAM: THICKNESS OF DROPLET AT APEX / RADC; R: NONDIMENSIONAL RADIAL POSITION 
PRINT 27 98 


BETA = 2: * RADO / (6. * V / PY) °° (1.7/3.3 
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PRINT *. ‘BETA: “7 SETA 

GAM = (4./ BETA*®*3 + (1. 
- (-4./ BETA*®*3 £€ 

DO 700 J = 1,NUMCOL 

Rie (25% Jager dt.) (2~ 

DEL(J) = ((1./GAM + GAM) 

MAKE DEL() DIMENSIONAL 

DEL(J) = RADO * DEL(J) 

CONTINUE 

RETURN 

END 


Chi ie Case F5 
ee.S) 


+ 16./ BETA**6) ** 5) 
(1. + 16./ BETA**6) oir: 
* NUMCOL) 
eye 


= 82) 289) =) Cle GAM 


SUBROUTINE UPGEOB (L, V, BO, VO, RO, TO, TR, NUMCOL, DEL, R, 
UPDATE THE HEIGHT OF THE DROPLET USING MODLE ‘B’ 

ALSO CALCULATE THE CURRENT CONTACT ANGLE THETA IN RADIANS 
THIS SUBROUTINE WRITTEN BY SUSAN TINKER 

AND TRANSLATED FROM QUICKBASIC TO FORTRAN BY GLENN WHITE 


fa) 


GAM) /2. 


THETA) 


DEFINE COMMON BLOCK TO ALLOW SUBROUTINE UPGEOB TO KEEP VALUES OF CERTAIN VARIABLES STATIC 


REAL VR, VL, SO, Z0, XO, AO, Li, L2, SR, AR, ZR, XR, RL, ZL, 
COMMON /GEOBBL/ VR, VL, SO, ZO, X0, AO, Li, L2, 

PARAMETER (PI = 3.14159265358979) 

REAL V, BO, VO, RO, TO, TR, DEL(12), R, THETA 


INTEGER NUMCOL 
CHARACTER*1 L 


SL 


SR, AR, ZR, XR, RL, ZL, SL 


eeeeeeeSECTION I - DETERMINATION OF INITIAL DROPLET PARAMETERS***rerere 


DEFINITION OF VARIABLES IN SECTION I 
TO: INITIAL CONTACT ANGLE; INPUT BY USER 
INITIAL DROPLET VOLUME; INPUT BY USER 
: D/D, WHERE D= DROPLET DIAMETER BEFORE IMPACT 
MAX DROPLET DIAMETER AFTER IMPACT, INPUT BY USER 
RECEDING ANGLE; INPUT BY USER 
: ANY DROP VOLUME BETWEEN 0 AND VO, 
: MAX DROPLET RADIUS AFTER IMPACT 
: DISTANCE FROM CENTER OF INITIAL ARC TO R-AXIS 
RADIUS OF DISK PART OF DROPLET 
: HEIGHT OF INITIAL DROPLET 
: DISTANCE ON R-AXIS FROM XO TO RPO 
AQ: DISTANCE ON R-AXIS FROM RPO TO R 
RPO: R-COORDINATE OF PO 
REAL A, C, CL, DO, DR, D, E, F, GAM, G, RPO, 
REAL SLP, S, T, Ul. UR, U2, VA1, VB1, VCl. VD1. 
REAL VA2, VB2, VC2, VD2, V2, X, XC, Yl, YR, Y2, 
INTEGER J 
IF (&% -EQ: °0°) GOTO 1500 
C = VO / (PI * RO ** 2) 
E = (3 * VO) / (PI © RO *e 2) 
SO «= (C + E) / 2. 
ZO = SO * {((1. / TAN(TO)) 
Yl = (ZO + SO) 
XO = RO - SO * ((1. 
Ul = (RO - XO) / Y¥1 
VAL = (PI * XO ** 2 © SO) 
VBI Ce) 2. eerie ((YL 225373) = CL phew oe Del 
vc1l = PI * Y1 ** 2 * XO © (ATAN(U1 / SQRT(1. - U1 ** 2)) 
(1. / 2.) © SIN(2. © (ATAN(U1 / SQRT(1. 
VD1 = PI * ZO * (RO f° 2 - XO ©* 2) 
Vil «© VAL + VB1 + VCl - VD 
IF ((V1l - VO) .GT. 18-12) 
E = SO 
GOTO 1010 
ELSEIF ((VO - V1) 
Cc = SO 
GOTO 1010 
ELSE 
GOTO 1020 


INPUT BY USER 


RPR, RP, R1, R2 
V1, VAR, VBR 
zc 


+ SQRT(1. © (1. / TAN(TO)) ** 2)) 


/ TAN(TO)) ¢ SQRT(1. + (1. / TAN(TO)) 


THEN 


-GT. 18-12) THEN 
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(COS(ATAN(U1 / SQRT(1. 


» RV, ROV 
, VCR, VDR 


© TAN(PI / 2. - TO) 


2)) 


= US s9)2)))9 


NOY Se? 2 diy ))) 


ee 3) 


END IF 

1020 AO = SO / TAN(TO) 
DO = SO * SQRT(1. + (1. / TAN(TO)) ** 2) 
RPO = RO - AO 


Gases eerverEND SECTION [ewereerrererersserereareraesecescoseeerersereresseT 
Cee: *eeeeerSECTION II - CALCULATE VL AND VReeeeeeeeoercerecereesererrerereTe 
C Reuetes TR: RECEDING ANGLE; INPUT BY USER 

(2b oe GE SR: HEIGHT OF DROPLET AT RECEDING SHAPE 

(SAE AR: DISTANCE ON R-AXIS FROM PR TO R 

Casaee RPR: R-COORDINATE OF POINT PR 

(Socio, cr VR: VOLUME OF DROPLET AT RECEDING ANGLE 

feterko. VL: VOLUME OF DROPLET AT LENS SHAPE 

yap RL: RADIUS AT LENS SHAPE IF NOT RO 

(5.0 cl SL: HEIGHT OF LENS 


L1 = RO - SO / TAN(TR) 
L2 = SQRT((SO / TAN(TR)) ** 2 + SO ** 2) 
IF (L2 .GE. L1) THEN 
SR = (RO / TAN(TR)) * ((SQRT(1. + (TAN(TR)) ** 2)) - 1.) 
AR = SR / TAN(TR) 
ZR = RO / TAN(TR) 
RPR = RO - AR 
VR la /G eee Pla aOR GMa a cee RO Mita oe GR eee) 
VL = VR 
SL = SR 
ZL = ZR 
ELSE 
SR = SO 
AR = SR / TAN(TR) 
RPR = RO - AR 
DR = SR * SQRT(1. + (1. / TAN(TR)) ** 2) 
ZR = (DR + AR) * TAN(PI / 2. - TR) 
XR = RO - DR - AR 
YR = ZR + SR 
UR = (RO - XR) / YR 
VAR. «= PI * XR ** 2 * SR 


VBR = 2. * PI * ((YR ** 3 / 3.) - (YR ** 3 / 3.) * (COS(ATAN(UR / SQRT(1. - UR ** 2)))) ©* 3) 
VCR = PI * YR ** 2 * XR * (ATAN(UR / SQRT(1. - UR ** 2)) 
es + (1. / 2.) © SIN(2. © (ATAN(UR / SQRT(1. - UR ** 2))))) 


VORB=) Pio * (ZReeGtRO ** 2) = XR ° S52) 
VR © VAR + VBR + VCR - VDR 
CL = RO - SO * ((1. / TAN(TR)) ¢ SQRT(1. + (1. / TAN(TR)) ** 2)) 
RL © RO - CL 
ZL = RL / TAN(TR) 
SL = SO 
VLeoe (Pies. 2) 2 1SR Ee (set RL 2? (202 eSR a2) 
C _@°eeeeEND SECTION [Jl eteeeeeceveseccccesereseseeseseseseeeeseeeereeeeerTseeTeerTHeneTeTeet 
C.... *°°*SECTION III - CONFIGURATION FOR ANY VOLUME BETWEEN 0 AND VO**%*eererecrcccrceserares 
1800 IF (V .EQ. VO) THEN 
(os eree INITIAL CONFIGURATION® eeeeceecccccovececcocensecororeseenteneseearereteHaeHaeTeTeTe 
DO 1001 J = 1, NUMCOL ‘ 
Xee RO © (2. * J = 22) / (22 Ss NONCOL) 
IF (X .LE. XO) THEN 
DEL(J) » SO 
ELSE 
DEL(J) = SORT((Z0 « SO) °° 2 - (X - XO) ** 2) - 20 
END IF 
2001 CONT INUE 
R = RO 
GAM = SO / (RO - XO) 
GOTO 1000 
END IF 
IF (V .GT. VL) THEN 
IF (V .GT. VR) THEN 
IF (L2 .GE. Li) THEN 
SLP = (SO - SR) / (AR - AO) 
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1030 


1040 


1002 


F = TR 

G = TO 

Teen (fs °G) 72; 

Als (AO eailten/ SLP)' © So) 7 (1. « (1. / SLB)O*)TAN(T))) 
S = A * TAN(T) 

D = SQRT(S ** 2+ A ** 2) 
xc = RO-D-A 

RP = RO-A 

zc = (D + A) / TAN(T) 

Y2 = (ZC + S) 

U2 = (RO - XC) / Y2 

VA2 = (PI * XC ** 2 © S) 


vB2 = 2. * PI * (Y¥2 ** 3 / 3.) * (1. - (COS(ATAN(U2 / SQRT(1. - U2 ** 2)))) ** 3) 


vC2 = PI * Y2 ** 2 * XC © (ATAN(U2 / SQRT(1. - U2 ** 2)) 
+ (1. / 2.) * SIN(2. * (ATAN(U2 / SQRT(1. - U2 ** 2))))) 
VD2 = PI * ZC * (RO ** 2 = XC ** 2) 
V2 = VA2 + VB2 + VC2 - VD2 
IF ((V2 - V) .GT. 1E-12) THEN 


G=T 
GOTO 1030 
ELSEIF ((V - V2) .GT. 1E-12) THEN 
F=T 
GOTO 1030 
ELSE 
GOTO 1050 
END IF 
ELSEIF (L1 .GT. L2) THEN 
F = TR 
G = TO 
Ti s2lFi eG) (<2 
S = SO 
A= S / TAN(T) 
De SORT(S ** 2¢A** 2) 


xC = RO-D-A 

RP = RO- A 

zC = (D + A) / TAN(T) 
Y2 = 20 +S 

U2 = (RO - XC) / Y2 

VAZ) =§PIs XC ** 2° S 


VB2 « 2. © PI * (Y2.** 3 / 3.) © (1. = (COS(ATAN(U2 / SQRT(1. - U2 ** 2)))) ** 3) 


VG2Ze=e Y2 92° 2 °F * XC * (ATAN(U2 / SOQRT(2. ~ U2 ** 2)) 
+ (1. / 2.) * SIN(2. * (ATAN(U2 / SQRT(1. - U2 ** 2))))) 
MOSee Lecce teORG 2 ree a KC en 2) 
V2 = VA2 + VB2 + VC2 - VD2 
EF UOMV2eeeV)) .GT.. LE-22) THEN 
GefT 
GOTO 1040 
ELSEIF ((V - V2) .GT. 1E-12) THEN 
F ef 
GOTO 1040 
ELSE 
GOTO 1050 
END IF 
ELSE 
CONTINUE 
END IF 
e@eereeeCONFIGURATION FOR VOLUME V > VReeeeeeecveceseceorerceoerereereneenersereeerETes 
DO 1002 J = 1, NUMCOL 
KepROeen(2. * J = 2%.) / C2). * NUMCOL) 
IF (X .LE. XC) THEN 
DEL(J) = S 
ELSE 
DEL (J) «= SQRT((ZC « &) 9° 2 = {X - KC) ©* 2) - ze 
END IF 
CONTINUE 
R = RO 
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SAM = S ;/ iRO - XC) 
ELSEIF (V .EQ. VR) THEN 
DO 1003 J = 1, NUMCOL 
X = RO ®* (2. * J - 1.) / (2. * NUMCOL) 
IF (X .LE. XR) THEN 
DEL(J) = SR 
ELSE 
DEL(J) = SQRT((ZR + SR) ** 2 - (X - XR) ** 2) - ZR 
END IF 
1003 CONTINUE 
R = RO 
GAM = SR / (RO - XR) 
ELSE 
Rl = RL 
R2 = RO 
1100 RV = (R1 + R2) / 2. 
S = SO 
A = S / TAN(TR) 
RP = RV-A 
D = S * SQRT(1. + (1. / TAN(TR)) ** 2) 
zc = (D + A) * TAN(PI / 2. - TR) 


xC = RV-D-A 
Y2 = ZC +S 


U2 = (RV - XC) / ¥2 
VAZee eri e* ecu" 2209S 
VB2 = 2. * PI * ((¥2 ** 3 / 3.) - (¥2 ** 3 / 3.) * (COS(ATAN(U2 / SQRT(1. - U2 ** 
vC2 = PI * ¥2 ** 2 * XC * (ATAN(U2 / SQRT(1. - U2 ** 2)) 
+ + (1. / 25)) © °SIN(2. * (ATAN(U2 / SQRT(1. - U2 ** 2))))) 
VD2 =EPL *gZC (* (RV e272 5- XC 7e* 2) 
V2 = VA2 + VB2 + VC2 - VD2 
IF ((V2 - V) .GT. 1E-12) THEN 
R2 = RV. 
GOTO 1100 
ELSEIF ((V - V2) .GT. 1E-12) THEN 
Rl = RV 
GOTO 1100 
ELSE 
GOTO 1200 
END IF 
Pe Ane @eerereCONFIGURATION FOR V < VR*teeeeercerreerererensreareerrerererereaeeteeeeneteTeteee 
1200 DO 1004 J = 1, NUMCOL 
X= RO * (2. *13 = 15) 0/7 (2. ) *2NOMCOL) 
IF (X .LE. XC) THEN 
DEL(J) = S 
ELSEIF (X .LE. RV) THEN 
DELAD) f= SORT ZCRoe SI it * F2- aX XC) 2882) 5 =" ZC 
ELSE 
DEL(J) = 0. 
END IF 
1004 CONTINUE 
R = RV 
GAM = S / (RV - XC) 
END IF 
ELSEIF (V .EQ. VL) THEN 
he ete eeeeeereCONFIGURATION FOR V @ VL® Pee ertrrrreencesarersecercesevereneresaeeTeree 
DO 1005 J = 1, NUMCOL 
Xe RO * (2. * J'= 1°) / (2. * NOMCOL) 
IF (X .LE. RL) THEN 
DEL(J) = SORT((ZL + SL) e* 2 - X ** 2) - ZL 
ELSE 
DEL(J) = O. 
END IF 
1005 CONTINUE 
R = RL 


GAM = SL / RL 


ELSE 
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2)))) 


** 3) 


_©eee8eeeCONFIGURATION FOR V < VLe eerececrreerecveeceesrrsegrrrrs 
Rov = (V / ((PI / 6.) * (1. - COS(TR)) © (3. © (SIN(TR)) ** 2 +929 -pCOS(TR)} ©" .2))) ** <2. 
RV = ROV * SIN(TR) 

S = ROV * (1. - COS(TR)) 
zC = ROV - S 
DO 1006 J = 1, NUMCOL 
X = RO * (2. © J - 1.) / (2. * NUMCOL) 
IF (X .LE. RV) THEN 
DEL(J) = SQRT((ZC + S) ** 2 - X ** 2) - ZC 
ELSE 
DEL(J) = 0. 
END IF 

1006 CONTINUE 
R = RV 
GAM = S / RV 

END IF 


eteneerEND SECTION [II**tererecceeseceecceverrrrereverrrererensaraenvererere 


1000 CONTINUE 
CALCULATE THETA FROM GAM, THE RATIO OF APEX TO DELTA(R) OF CURVED SURFACE 


THETA = ATAN(((1./GAM + GAM) **2 /4. - 1.) ** -.5) 
RETURN 
END 


SUBROUTINE WEIGHT (NUMWTI, NUMNOD, OT, RADO, W) 
CALCULATE THE WEIGHT MATRIX W CAREFULLY FOR THE FIRST NUMWTI TIME STEPS 
WARNING: GOOD ONLY FOR DTSHORT = 0.1 AND DTLONG 2 1.0 !!!!!!! 
PARAMETER (PI = 3.14159265358979) 
INTEGER NUMWTI, NUMNOD 
REAL DT, RADO, W(78, 78, 10) 
Pea can NUMWTI: NUMBER OF TIME STEPS TO BE PRECALCULATED 
REAL ALPHA, DTO, TO, CONST1, DRO, RO, R. ARG1, ARG2, ARG3, SUM 
INTEGER I, J, K, L, LMAX 
Aa Ss ALPHA: THERMAL DIFFUSIVITY OF SOLID; DTO: RECOLLECTION TIME INTERVAL; 
Amerie TO: RECOLLECTION TIME; CONST1: ALGEBRAIC EXPRESSION FROM GREEN'S FUNCTION; 
einer DRO: DUMMY RADIAL POSITION INTERVAL, RO: DUMMY RADIAL POSITION; 
aopeet te R: RADIAL POSITION; ARG1. ARG2, ARG} GREEN'S FUNCTION ARGUMENTS; 
suena ene SUM: INTEGRATION SUMMATION; L: INTEGRATION INDEX: LMAX: NUMBER OF STEPS PER TIME INTEGRATION 
REAL ERF, BESIOE 
1 eee ERF: ERROR FUNCTION; BESIOE: EXPONENTIAL BESSEL FUNCTION IO(ARG) X EXP(-ARG) 
sites 7% (FROM LINKED LIBRARY) 
IF (DT.EQ.0.1) LMAX e@ 25 
IF (DT.EQ.1.0) LMAX e 250 
ALPHA = 5.79E-07 
OTO = OT 
CONST1 © (4 * PI * ALPHA) °° - § 
DO 750 K e 1,NUMWTI 
DO 760 I = 1,78 
pO 770 J 2 1.78 
SUM = OQ. 
IF (I LS. 66) THEW 
DRO «= RADO / 12 
RO e RADO © (I - $) / 143 
ELSEIF (I LE 66) THD 
DRO © RADO / 6 
RO © RADO © (@ © (1:2 - 686 - § 6 
ELSEIF ‘(I LE 7S) TH 
DRO »« RADO / 3} 


RO © RADO © (7 ¢ 1:2 - 66: - § ,) 
ELSE 

DRO e RADO 

RO © RMADO ® (10 © (if Se Wo | 
DOD IF 


IF (J .LE. 46) THEW 
Re RADO ® (J - $) 32 
ELSEIF (J LE 66) THE 
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R = RADO © (4. + ((J - 48) - .5i ¢ §.3 
ELSEIF (J .LE. 75) THEN 


R = RADO * (7. + ((J - 66) - .5) / 3.) 
ELSE 

R= RADO * (10. + ((3 > 75) =" -5)) 
END IF 


IF (I .EQ. J) THEN 

DO 780 L = 1,LMAX 

To = (K-1+tUL/ (1. * LMAX)) * DTO 

ARG1 = R * RO / (2. * ALPHA * TO) 

ARG2 = (DRO / 2.) / (4. * ALPHA * TO) ** .5 

780 SUM = SUM + TO ** -1 * BESIQE(ARG1) * ERF(ARG2) 
W(J, I, K) = R * SUM * DTO / (1. * LMAX) 
PRINT *, J, I, K, WJ, I, K) 

ELSE 

DO 790 L = 1, LMAX 

TO 0, (K = 01 7b" 7/7 (12) ®* LMAX))*2°sDTO 
ARG1 = R * RO / (2. * ALPHA * TO) 
ARG3 = -(R - RO) ** 2 / (4. * ALPHA * TO) 


790 SUM = SUM + TO ** -1.5 * BESIOE(ARG1) * EXP(ARG3) 
W(J, I, K) = CONST1 * RO * DRO * SUM * DTO / (1. * LMAX) 
END IF 
770 CONTINUE 


760 CONTINUE 
750 CONTINUE 
RETURN 
END 


REAL FUNCTION WGHT(J, I, K, DT, RADO, ®) 
nado FUNCTION TO GIVE THE PROPER WEIGHT FUNCTION VALUE 
0 as WARNING: GOOD ONLY FOR DTSHORT = 0.1 AND DTLONG = 1.0 !!!!!!! 
PARAMETER (PI = 3.14159265358979) 
INTEGER J, I, K 
REAL DT, RADO, W(78, 78, 10) 
REAL ALPHA, DTO, TO, CONST1, DRO, RO. R, ARG1, ARG2, ARG3 
INTEGER NUMWTI 
eaeiate ALPHA: THERMAL DIFFUSIVITY OF SOLID; DT0: RECOLLECTION TIME INTERVAL; 
Sats TO: RECOLLECTION TIME; CONST1: ALGEBRAIC EXPRESSION FROM GREEN’S FUNCTION; 
Aeeyercers DRO: DUMMY RADIAL POSITION INTERVAL; RO: DUMMY RADIAL POSITION; 
Berle R: RADIAL POSITION; ARG1, ARG2, ARG3: GREEN’S FUNCTION ARGUMENTS; 
Societe NUMWTI: NUMBER OF TIME STEPS TO BE PRECALCULATED 
REAL ERF, BESIOE 
Sieteters ERF: ERROR FUNCTION; BESIOE: EXPONENTIAL BESSEL FUNCTION IO(ARG) X EXP(-ARG) 
Sorts 7 (FROM LINKED LIBRARY) 
IF (DT.EQ.0.1) NUMWTI = 10 
IF (DT.EQ.1.0) NUMWTI = 2 
IF (K.LE.NUMWTI) THEN 
WGHT = W(J, I, K) 
RETURN 
END IF 
ALPHA = 5.79E-07 
DOTO = OT 
CONST1 = (4 * PI * ALPHA) °° -.5 
TOR@e(Ke= © 5) *ODTO 
IF (I .LE. 48) THEN 
DRO = RADO / 12. 
RO = RADO * (I - .S) / 12. 
ELSEIF (I .LE. 66) THEN 
DRO = RADO / 6. 
RO = RADO * (4. + ((I - 48) - $) / 6) 
ELSEIF (I .LE. 75) THEN 
DRO = RADO / 3. 


ROGeDRADO © (7o 0s (CE 6-766) == Si) 795) 
ELSE 
DRO = RADO 


A22 


RO'’s RADO © :10. <= {( 
END IF 
IF (J .LE. 48) THEN 
Ros RADO *9 (5 = .8) / 22. 
ELSEIF (J .LE. 66) THEN 
R = RADO * (4. + ((J - 48) - .5) / 6.) 
ELSEIF (J .LE. 75) THEN 


HH 
‘ 
4 
uw 
‘ 
ut 


Res RADO) 2 (75 +) (ld (=) 66) 5-25) 0/7 3s) 
ELSE 

R= RADO "5 (10. +1 =" 75) = -$)) 
END IF 


ARG1 = R * RO / (2. * ALPHA * TO) 
IF (I .EQ. J) THEN 
ARG2 = (DRO / 2.) / (4. * ALPHA * TO) ** .5 
WGHT = R * TO ** -1 * BESIOE(ARG1) * ERF(ARG2) * DTO 
ELSE 
ARG3 = -(R - RO) ** 2 / (4. * ALPHA © TO) 
WGHT = CONST1 * RO * TO ** -1.5 * BESIOE(ARG1) * EXP(ARG3) * DRO * DTO 
END IF 
RETURN 
END 


REAL FUNCTION XI (TI, PSAT) 
FUNCTION TO GIVE XI GIVEN TEMPERATURE 
XI: MOLAR FRACTION OF WATER VAPOR 
REAL TI, PSAT(0:100) 
INTEGER TSAT 
REAL TVAL, XA, XB 
TSAT = -1 
TVAL = TI 
IF (TVAL .GT. 100.) TVAL = 100. 
900 CONTINUE 
TSAT «= TSAT + 1 
IF (TSAT .LT. TVAL - 1.) GOTO 900 
XA = PSAT(TSAT) / 1.01325 
XB = PSAT(TSAT + 1) / 1.01325 
yet ate LINEARLY INTERPOLATE XI 
XI = XA + (XB - XA) * (TI - TSAT) 
RETURN 
END 
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Appendix C Material Properties as a Function of Temperature 
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FIGURE 48 
Lewis Number of Steam in Air vs. Temperature 
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FIGURE 49 


Constant Pressure Specific Heat for Air Cp,air VS. Temperature 
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FIGURE 50 
Molar Fraction of Water Vapor at Liquid-Vapor Interface xj vs. Temperature 
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FIGURE 51 
Ratio of A/ Vfg for Water vs. Temperature 
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FIGURE 52 


Thermal Conductivity of Water k; vs. Temperature 
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FIGURE 53 


Thermal Diffusivity of Water a vs. Temperature 
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y = 999.79 + 6.4738e-2x - 8.2228e-3x*2 + 5.8548e-5x*3 - 3.2774e-7x*4 + 8.5038e-10x45 
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FIGURE 54 


Density of Water p; vs. Temperature 
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TABLE 2 
Thermodynamic Properties of Saturated Water [23] 


t > U ) 10°3 vg g g g/VIg 

Fe] wea || emeangy| cages] crv) | ims) | Ga | tnd) 
; TO) OUU 100 | 0.0000 te 1.0002 | 999.800 | 2062/5 | U0. 77/' J} 2501.4 | 12.1264 | 
—710-0006567 [0.0064811[ 1.0002 | 998.800 | 192577 [192.5760] 2499.0 | 12.9767, 
2 [0.0007056 [0.0069637| 1.0001 | 999.900 [179889 |179.8880| 2496.7 | 13.8792 
T 3 [0.007577 |0.0074779| 1.0001 | 999.900 | _ 168132 |168.1310} 2494.3 | 14.8355 | 
-—410,0008131 100080247] 1.0001 | 999.900 |_157232 | 157.2310| 2491.9] 15.8487 | 
5 10,0008721 [0.086070] 1.0001 | 999.900 | 147120 [147.1190] 2489.6 | 16.9224 
6 [0.0009349 | 0.0082267] 1.0001 | 999.900 |" 137734 [137.7330 | 2487.2 | 16.0581 
T 8 [0.0010724|0.0105838] 1.0002 | 999.800 | 120917 |120.9160| 2482.5 | 20.5308 | 
T 9 [0.0011477[0.0113269| 1.0003 | 999.700 |_113386 |113.3850| 2480.1 | 21.8733 | 
| t2 [0.0014022]0.0138386] 1.0005 | 999.500 | 93784 | 93.7830 | 2473.0 | 26.3694 | 
| 15 [0.0017051 | 0.0168280| 1.0009 | 999.101 | 77926 | 77.9250 | 2465.9 | 31.6445 | 
[18 [0.0020640[0.0203701| 1.0014 | 998.602 | 65038 | 65.0370 | 2458.8 | 37.6062, 
| 20] 0.002339 | 0.023084 |__ 1.0018 | 998.203 | 57791 | 57.7900 | 2454.1 | 42.4658 | 
| 24] 0.002985 | 0.029460] 1.0027 | 997.307 | 45883 | 45.8820 | 2444.7 | 53.2823 | 
| 25 | 0.003169 | 0.031276 | 1.0029 | 997.108 | 43360 | 43.3590 |_ 2442.3 | 56.3274 | 
| 26 | 0.003363 | 0.033190| 1.0032 | 996.810 | 40994 | 40.9930 | 2439.9 | 59.5199 | 
28 | 0.003782 | 0.037325 | 1.0037 | 996.314 | 36690 | 36.6800 | 2435.2 | 66.3741 | 
| 29 | 0.004008 | 0.039556 | 1.0040 | 996.016 | 34733 | 34.7320 | 2432.8 | 70.0449 | 
| 30 | 0.004246 | 0.041905] 1.0043 | 995.718 | 32894 | “32.8930 | 2430.5 | 73.8911 | 
| 31 | 0.004496 | 0.044372] 1.0046 | 995.421 | 31165 | 31.1640 | 2428.1 | 77.9136 | 
| 32 | 0.004759 | 0.046968] 1.0050 | 995.025 | 29540 | 29.5390 | _ 2425.7 | 82.1186 | 
| 33] 0.005034 | 0.049682 | 1.0053 | 994.728 | 28011 | 28.0100 | 2423.4 | 86.5191 | 
| 34] 0.005324 | 0.052544 | 1.0056 | 994.431 | 26571 | 26.5700 | 2421.0 | 91.1178 | 
| 36 | 0.005947 | 0.058692] 1.0063 | 993.739 | 23940 | 23.9390] 2416.2 | 100.932 | 
| 38 | 0.006632] 0.065453] 1.0071 | 992.950 | 21602 | 21.6010 | 2411.5 | 111.638 | 
| 40 | 0.007384 | 0.072874] 1.0078 | 992.260] 19523 | 19.5220] 2406.7 | 123.281 | 


41 | 0.007786 | 0.076842 991.867 18.5690 
42 0.008208 1.0086 | 991.473 17671 | 17.6700} 2401.9] 135.931 
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TABLE 2 (Continued) 
Thermodynamic Properties of Saturated Water [23] 
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-o.020881 | 0.206079] 1.0177 | 962.608 | 7346 7.3450 | 236.0 | 320.763. 
62 
63 
65] 0.02503] 0.24703] 1.0199 | 980.488 | 6197] 6.1960 | 2346.2 | 376.665 | 
66} 0.02617 | 0.25828 | 1.0205 | 979.912] 5943] 5.9420 | 2343.7] 304.431 | 
68] 0.02859 0.28216 | 1.0217 | 978.761| 5471] 5.4700| 2398.8 | 427.570 
69] 0.02986 | 0.20470 | 1.0222 | 978.282 | 5252] 5.2510 | 2336.3 | 444.927 | 
71] 0.03256 | 0.32134 [ 1.0234 | 977.135] 4843] 4.8420 | 2331.4 | 481.408 | 
7] 0.04192 | 0.41372] 1.0272 | 973.520 3822] 3.8210 | 2316.3 | 606.207 | 
(79 [0.04550 | 0.44905] 1.0285 | 972.200 | 3539] 3.5380] 2311.3 | 653.204 | 
(a0 | 0.04739] 0.46770 1.0291 | 971.723 | 3407] 3.4060 | 2308.8 | 677.868 | 
82] 0.05136 | 0.50688 | 1.0305 | 970.403] 3160] 3.1590] 2303.7] 729.257) 
83] 0.05345] 0.52751 | 1.0311 | 969.838 | 3044] 3.0430 | 2301.1 | 756.202 | 
84] 0.05560 | 0.54873 1.0318 | 969.180 | 2934] 2.9330] 22986 | 783711 | 
85; 0.05783} 0.57074| 1.0325 | 968.523 2296.0 | 812.178 
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Thermodynamic Properties of Saturated Water [23] 


TABLE 2 (Continued) 
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TABLE 3 
Thermal Properties of Saturated Water [16] 
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Appendix D.1 Output Data: Effect of Shape Factor B 
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FIGURE 55 
Transient Droplet Volume: Runs 1-9 (Effect of Bo) 
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FIGURE 56 HIE 
Transient Rate of Change of Droplet Volume: Runs 1-9 (Effect of Bo) 
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Transient Contact Angle @: Runs 1-9 (Effect of Bo) 
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FIGURE 58 TIME (s) 
Transient Shape Factor B: Runs 1-9 (Effect of Bo) 
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BOTTOM SURFACE AVERAGED FLUX (kW/m*2) 
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FIGURE 59 Me 
Transient Area Averaged Solid-Liquid Flux: Runs 1-9 (Effect of Bo) 
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DROPLET SURFACE AVERAGED TEMPS. (C) 
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FIGURE 60 

Transient Upper and Lower Surface Averaged Droplet Temperatures: Run 1 
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DROPLET SURFACE AVERAGED TEMPS. (C) 
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FIGURE 61 TIME (s) 
Transient Upper and Lower Surface Averaged Droplet Temperatures: Run 2 
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DROPLET SURFACE AVERAGED TEMPS. (C) 


O; O08 5; @b10) ob15; 2620) ye2d: ac, 9e35. 2740, 9:45. 2 90. 255° «60 
TIME (s) 

FIGURE 62 

Transient Upper and Lower Surface Averaged Droplet Temperatures: Run 3 
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DROPLET SURFACE AVERAGED TEMPS. (C) 
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FIGURE 63 | TIME (s) 
Transient Upper and Lower Surface Averaged Droplet Temperatures: Run 4 
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DROPLET SURFACE AVERAGED TEMPS. (C) 
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FIGURE 64 TIME (s) 

Transient Upper and Lower Surface Averaged Droplet Temperatures: Run 5 
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DROPLET SURFACE AVERAGED TEMPS. (C) 
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FIGURE 65 TIME (s) 
Transient Upper and Lower Surface Averaged Droplet Temperatures: Run 6 
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DROPLET SURFACE AVERAGED TEMPS. (C) 
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FIGURE 66 eS) | 
Transient Upper and Lower Surface Averaged Droplet Temperatures: Run 7 
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DROPLET SURFACE AVERAGED TEMPS. (C) 
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FIGURE 67 des) 
Transient Upper and Lower Surface Averaged Droplet Temperatures: Run 8 
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DROPLET SURFACE AVERAGED TEMPS. (C) 
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FIGURE 68 UME 6) 
Transient Upper and Lower Surface Averaged Droplet Temperatures: Run 9 
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RECOLLECTION MEMORY (s) 
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TIME (s) 


FIGURE 69 
Transient Recollection Memory: Runs 1-9 (Effect of Bo) 
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Appendix D.2 Output Data: Effect of Initial Volume Vj 
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FIGURE 70 TIME (s) 
Transient Droplet Volume: Runs 5, 10, 11 (Effect of Vo) 
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FIGURE 71 Me 8) 
Transient Rate of Change of Droplet Volume: Runs 5, 10, 11 (Effect of Vo) 
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FIGURE 72 TIME (s) 
Transient Contact Angle 8: Runs 5, 10, 11 (Effect of Vo) 
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FIGURE 73 TIME (s) 
Transient Shape Factor B: Runs 5, 10, 11 (Effect of Vo) 
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BOTTOM SURFACE AVERAGED FLUX (kW/m/2) 
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FIGURE 74 

Transient Area Averaged Solid-Liquid Flux: Runs 5, 10, 11 (Effect of Vo) 
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DROPLET SURFACE AVERAGED TEMPS. (C) 
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FIGURE 75 TIME (s) 
Transient Upper and Lower Surface Averaged Droplet Temperatures: Run 10 
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DROPLET SURFACE AVERAGED TEMPS. (C) 
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FIGURE 76 TIME (s) 
Transient Upper and Lower Surface Averaged Droplet Temperatures: Run 11 
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Appendix D.3 Output Data: Effect of Initial Temperatures 
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Transient Droplet Volume: a) Run 15; b) Run 12; c) Run 13; d) Run 5; 
e) Run 14 (Effect of Initial Temperatures) 
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FIGURE 78 TIME (s) 
Transient Rate of Change of Droplet Volume: a) Run 15; b) Run 12; c) Run 13; 


d) Run 5; e) Run 14 (Effect of Initial Temperatures) 
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CONTACT ANGLE THETA (deg) 
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FIGURE 79 TIME (s) 

Transient Contact Angle @: a) Run 15; b) Run 12; c) Run 13; d) Run 5; 

e) Run 14 (Effect of Initial Temperatures) 
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FIGURE 80 TIME (s) 
Transient Shape Factor B: a) Run 15; b) Run 12; c) Run 13; d) Run 5; e) Run 14 
(Effect of Initial Temperatures) 
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BOTTOM SURFACE AVERAGED FLUX (kW/m‘2) 


25 
TIME (s) 


Transient Area Averaged Solid-Liquid Flux: a) Run 15; b) Run 12; c) Run 13; 
d) Run 5; e) Run 14 (Effect of Initial Temperatures) 


FIGURE 81 


150 


Appendix D.4 Output Data: Effect of Geometry Model 
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FIGURE 82 TIME (s) 
Transient Droplet Volume: d) Run 5; f) Run 17; g) Run 20; h) Run 18; 
i) Run 16; j) Run 19 (Effect of Droplet Geometry) 


0.60 


wee ge wwccobogecccccccboscc ccc cc chores ccc cet seccsccecccbsceccccsesbesccccccccbocoescoecoecbeccccceccebosccocess= 
. 


weccce 


ee eee 
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FIGURE 83 TIME (s) 

Transient Rate of Change of Droplet Volume: d) Run 5; f) Run 17; g) Run 20; 
h) Run 18; i) Run 16; j) Run 19 (Effect of Droplet Geometry) 
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FIGURE 84 TIME (s) : 
Transient Contact Angle @: d) Run 5; f) Run 17; h) Run 18; i) Run 16; 

j)) Run 19 eee of ela eee 
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FIGURE 85 NONDIMENSIONAL TIME (ttevap) 
Contact Angle @ vs. Time Scaled by t: d) Run 5; f) Run 17; h) Run 18; 
i) Run 16; j) Run 19 (Effect of Droplet Geometry) 
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FIGURE 86 TIME (s) 
Transient Shape Factor B: d) Run 5; f) Run 17; g) Run 20; h) Run 18; i) Run 16; }) 
Run 19 (Effect of Droplet Geometry) 
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FIGURE 87 TIME (s) 


Transient Area Averaged Solid-Liquid Flux: d) Run 5; f) Run 17; g) Run 20; 
h) Run 18; i) Run 16;7j) Run 19 (Effect of Droplet Geometry) 
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Appendix E FORTRAN Code for Constant Heat Flux Model 


ten, 


= SLENN WHITE 

Cc (301) 405-5334 

Cc ADVISOR: MARINO DIMARZO (301) 405-5257 

S MASTER'S THESIS -- UNIVERSITY OF MARYLAND AT COLLEGE PARK 

(S EVAPORATIVE COOLING WITH RADIANT HEAT INPUT 

c WRITTEN FOR NIST CENTER FOR FIRE RESEARCH 

a 

c SUDDEN PLACEMENT AND REMOVAL OF DISK HEAT SINK ON A SEMI-INFINITE SOLID 
S 

C234567 


PROGRAM DISK 
REAL RADO, TIME(24), TEVAP, T, TSO, KS, ALPHA, R, QD, QI, ARG1, ARG2, ARG3 


INTEGER I, J 
(Se tn RADO: RADIUS OF WETTED AREA; TIME: TIME SINCE INTRODUCTION OF DISK 
TEVAP: EVAPORATION TIME; T: SURAFCE TEMPERATURE; TSO: INITIAL SURFACE TEMPERATURE; 
KS: CONDUCTIVITY OF SOLID; ALPHA: THERMAL DIFFUSIVITY OF SOLID; 
Cert R: RADIAL POSITION; QD: DISK FLUX; QI: INITIAL FLUX; 
Sae6.0 G ARG1, ARG2, ARG3: ARGUMENTS TO SOLUTION INTEGRAL 


Cian = SPECIFY VARABLES TO BE USED BY SUB QAGI 
REAL EPSABS, EPSREL, RESULT, ABSERR, BOUND, WORK(5000), F 
INTEGER INF, IER, LIMIT, LENW, IWORK(1000), NEVAL, LAST 
EXTERNAL F 
COMMON ARG1, ARG2, ARG3 
OPEN (30, FILE = ‘valdiski.dat’) 


Cyaere a DEFINE PHYSICAL PARAMETERS 
RADO = 2.673009E-3 
TSO = 130. 
QD = 36534. 
QI = 2578. 
KS = 1.297 
ALPHA = 5.79E-7 
TEVAP = 31.5 


fests. ghekcins DEFINE THE TIMES TO BE USED 
DATAST IME) | /0.27, 0 iz teenie syeeeii ed ce Se Snel Oe Slane 2 One 2 Sees Oo 
° Jepsen Sar Fae, S957 40n 5 450, 501, 55260.) 


Creme DEFINE VARIABLES NEEDED BY SUB QAGI (CMLIB LIBRARY) USED FOR CALCULATION OF SEMI-INFINITE INTEGRAL 


BOUND = .00001 
INF = 1 

EPSABS = .001 
EPSREL = .001 
LIMIT = 1000 
LENW = 5000 


(Sp6 hos DO LOOP FOR CALCULATING THE DIMENSIONLESS TEMPERATURE AS A FUNCTION OF R 
do 15 I= 1, 24 
Do 10 J = 1, 78 
IF (J .LE. 48) THEN 
RessRADOT@E (0S <9 05) 97) Lae 
ELSEIF (J .LE. 66) THEN 


R = RADO * (4. + ((J - 48) - .5) / 6.) 
ELSEIF (J .LE. 75) THEN 

Rw PRADOP S27 ee eC (ure 6G ee S30 3'e) 
ELSE 

R = RADO * (10. + ((J - 75) - .§)) 
END IF 


ARG1 = R / RADO 
ARG2 = (ALPHA * TIME(I)) ** .5 / RADO 
IF (TIME(I).LE.TEVAP) THEN 
ARG3 = 0. 
ELSE 
ARG3 = (ALPHA * (TIME(I) - TEVAP)) °* .5 / RADO 
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CALL QAGI ‘F, BOUND, INF, EPSABS, EPSREL, RESULT, ABSERR,NEVAL, IER, LIMIT, LENW, LAST, WORK, WORK) 
T = TSO - (QD + QI) * RADO * RESULT / KS 
PRINT *, TIME(I), J, T 
WRITE (30, S) TIME(I), J, R / RADO, T 
10 CONTINUE 
1S CONTINUE 
5 FORMAT (F8.5, IS, 2F15.6) 
STOP 
END 


REAL FUNCTION F(X) 

REAL X 

REAL ARG1, ARG2, ARG3 

REAL ERF, BESIOE 

ERF: ERROR FUNCTION; BESJO, BESJ1: BESSEL FUNCTIONS JO AND Jl 

(FROM LINKED LIBRARY) 

COMMON ARG1, ARG2, ARG3 

F = BESJO(ARG1 * X) * BESJ1(X) © (ERF(ARG2 * X) - ERF(ARG3 * X)) / X 
RETURN 

END 
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